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Chapter  1 
Introduction 


1.1  Review  of  Phase  I  Achievements 

Significant  milestones  have  been  reached  through  the  work  of  Phase  I.  These  are  in  four 
related  categories.  Firstly,  the  modeling  and  prediction  of  flutter  and  LCO  of  an  airfoil 
and  control  surface  with  structural  freeplay  in  2D  transonic  flow  has  been  accomplished. 
Secondly,  the  modeling  and  prediction  of  flutter  and  LCO  of  an  airfoil  undergoing  single 
degree  of  freedom  pitch  motions  due  to  nonlinear  aerodynamic  effects  associated  with  large 
(inviscid)  shock  motions  in  2D  flow  has  been  achieved.  Thirdly,  the  modeling  and  predic¬ 
tion  of  flutter  of  a  3D  flow  about  an  elastic  wing  has  been  accomplished.  Note  that  in  the 
latter  category,  we  are  modeling  the  nonlinear  static  or  steady  flow  effects  in  the  transonic 
regime  even  though  the  aerodynamic  model  is  dynamically  linear,  i.e.,  the  shock  motions  are 
assumed  to  be  proportional  to  the  airfoil  motion  in  this  model.  And  finally,  a  2D  numerical 
simulation  of  the  NLR  7301  supercritical  airfoil  limit  cycle  oscillations  was  performed  using 
the  time-marching  procedure  of  the  CFL3D  code.  This  fourth  category  of  work  is  to  estab¬ 
lish  benchmarks  for  the  viscous  POD/ROM  methodology  developed  in  Phase  I  and  to  be 
developed  in  Phase  II. 


1.2  Phase  II  Goals  and  Technical  Objectives 

The  overarching  goal  for  the  Phase  II  work  is  to  build  on  the  very  substantial  progress 
made  in  Phase  I  and  to  develop  a  capability  for  a  highly  computationally  efficient  and  phys¬ 
ically  accurate  mathematical  modeling  of  limit  cycle  oscillations  (LCO)  and  other  nonlinear 
aeroelastic  phenomena.  The  approach  uses  the  concepts  of  aerodynamic  modes  as  well  as 
structural  modes.  Such  models  provide  substantially  improved  physical  understanding  and 
also  more  accurate  prediction  of  LCO  in  particular  and  nonlinear  aeroelastic  responses  in 
general.  This  capability  may  also  lead  to  a  rational  prediction  of  buffet  onset  due  to  aero- 
dynamically  nonlinear  oscillations. 

Phase  II  technical  objectives  included: 

•  Develop  the  2-D  N-S  version  of  the  above  proposed  approach 

•  Investigate  further  the  initial  condition  impart  on  LCO  via  a  2-D  time-domain  method 
(CFL3D/N-S  code) 
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•  Extend  the  geometry  capability  to  3-D  including  wing/store  configurations 

•  Generalize  the  2-D/3-D  methods  to  include  Harmonic  Balance  technique  for  nonlinear 
LCO  analysis 

•  Include  multiple  structural  mode  analysis  capability 

•  Perform  LCO  case  studies  including  wing-alone  and  F-16  wing/store  systems. 

1.3  Proposed  Method  for  Phases  I/II:  A  Brief  Review 

1.3.1  Limit  Cycle  Oscillations 

Limit  Cycle  Oscillations  (LCO)  are  known  to  occur  on  various  operational  vehicles.  In  par¬ 
ticular,  LCO  has  been  a  prevalent  aeroelastic  problem  on  several  current  fighter  aircraft.  It 
usually  occurs  for  aircraft  with  external  stores  throughout  the  transonic  flight  regime  (Refs 
[1,  2,  3,  4]).  Complicated  by  the  aircraft-store  system,  its  mechanism  still  remains  to  be  fully 
understood.  Meanwhile,  a  business  jet  wing  LCO  was  also  reported  recently  (Refs  [5,  6]). 
This  has  been  a  serious  concern  since  there  are  few  analytical  techniques  available  for  LCO 
prediction  and  an  insufficient  understanding  of  its  physics. 

LCO  can  be  characterized  as  sustained  periodic  oscillations  which  neither  increase  or  de¬ 
crease  in  amplitude  over  time  for  a  given  flight  condition.  Using  an  s-domain  unsteady 
aerodynamic  model  of  the  aircraft  and  stores,  Chen,  Sarhaddi  and  Liu  (Ref  [7])  have  shown 
that  wing/store  LCO  can  be  a  post-flutter  phenomenon  whenever  the  flutter  mode  contains 
low  unstable  damping.  This  type  of  flutter  mode  is  called  a  ’’hump  mode” .  Since  the  aircraft 
structure  usually  contains  structural  nonlinearity,  such  as  friction  damping,  this  amplitude- 
dependent  friction  damping  can  suppress  the  growth  of  amplitude,  thus  resulting  in  a  steady 
state  oscillation.  This  is  known  as  the  nonlinear  structural  damping  (NSD)  model  of  the 
wing/store  LCO.  Although  not  thoroughly  proven  through  tests  or  numerical  simulations, 
results  of  the  NSD  show  excellent  correlation  with  flight  test  LCO  data  of  F-16  throughout 
subsonic  and  transonic  Mach  numbers. 

On  the  other  hand,  other  researchers,  notably  Cunningham  and  Meijer  (Ref  [8])  believe 
that  the  wing/store  LCO  is  due  largely  to  transonic  shock  oscillation  and  shock  induced 
flow  separation.  This  is  called  the  Transonic  Shock/Separation  (TSS)  model.  Edwards  has 
suggested  the  TSS  model  and  viscous  effects  are  two  major  factors  that  cause  transonic  LCO 
for  wings.  He  also  has  studied  the  shock  buffet  phenomenon  in  addition  to  transonic  LCO 
(Ref  [9]).  It  should  be  noted,  however,  that  there  is  no  conflict  between  the  NSD  model  and 
the  TSS  model  in  that  both  physical  effects  may  contribute  to  LCO. 

Note  that  the  method  in  Ref  [7]  used  a  damping  correlation  technique  for  LCO-onset  predic¬ 
tion,  whereas  that  of  Ref  [8]  is  a  semi-empirical  method  that  requires  steady  and  unsteady 
data  input  from  wind  tunnel  measurement. 
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1.3.2  LCO  Prediction  Methods 


Recent  renewed  interest  in  LCO  is  motivated  by  the  need  to  better  predict  and  understand 
fighter  LCO  and  also  by  the  rapid  advent  of  CFD  methodology  in  aeroelasticity.  Cur¬ 
rently,  there  are  two  potential  computational  approaches  for  LCO  prediction/investigation: 
the  time-marching  approach  using  a  high-level  CFD  code  such  as  CFL3D  (Ref  [10])  devel¬ 
oped  and  supported  by  NASA/Langley,  and  the  Frequency-Domain  POD/ROM  EigenMode 
approach  (Ref  [11])  originated  by  Dowell  and  Hall  of  Duke  University.  The  former  is  a  con¬ 
ventional  time-domain  CFD  method  whereas  the  latter  is  a  frequency-domain  CFD  method 
using  aerodynamic  eigenmodes. 

For  nearly  twenty  years,  the  aerospace  industry  has  been  lacking  of  a  viable,  efficient  CFD 
method  for  transonic  flutter  applications.  The  time-marching  CFD  unfortunately  remain 
inefficient  as  a  practical  tool.  Under  an  AFOSR/STTR  Phase  I  contract,  ZONA  has  re¬ 
cently  conducted  a  feasibility  study  on  the  simulation  of  transonic  LCO  for  a  supercritical 
NLR7301  airfoil  using  the  CFL3D/Navier-Stokes  code.  It  was  found  that  the  time-marching 
LCO  solution  is  turbulence-model  dependent  and  time-step  sensitive.  One  LCO  solution 
case  typically  would  take  96  hours  on  a  1GHz  CPU  computer. 

Earlier  work  of  Bendiksen  (Ref  [12])  and  recent  work  of  Sheta  et  al  (Ref  [13])  also  used 
time-marching  CFD  approach  for  Euler  transonic  and  N-S  incompressible  LCO  simulations, 
respectively.  Based  on  the  cases  studied  by  ZONA,  the  time-marching  Euler  simulation  of 
a  transonic  LCO  requires  one  half  of  the  computing  time  required  of  a  Navier-Stokes  simu¬ 
lation,  due  to  the  reduction  in  grid  points.  The  computing  time  required  for  incompressible 
Navier-Stokes  LCO  simulation  is  about  two  orders  of  magnitude  less  than  a  transonic  Navier- 
Stokes  LCO  simulation,  due  to  the  absence  of  transonic  nonlinearity,  thus  allowing  much 
coarser  time  steps  and  a  reduction  in  subiterations.  For  transonic  Navier-Stokes  simulation 
of  LCO,  the  requirements  in  time  steps,  grid,  subiterations,  etc.,  become  much  more  strin¬ 
gent,  resulting  in  days  of  computing  time  for  one  LCO  case.  This  leads  to  the  conclusion 
that,  if  a  transonic  Navier-Stokes  simulation  were  demanded,  a  time-marching  approach 
would  be  computationally  inefficient  as  a  viable  approach  for  realistic  3D  transonic  LCO 
prediction/investigations.  On  the  other  hand,  the  Frequency-Domain  POD/ROM  (Proper 
Orthogonal  Decomposition  /  Reduced  Order  Model)  method  of  Duke  University  is  consid¬ 
ered  a  recent  breakthrough  in  the  transonic  computational  aeroelasticity.  Within  the  last  five 
years,  the  Frequency-Domain  POD/ROM  method  of  Dowell  and  Hall  has  been  proven  to  be 
highly  efficient  and  shown  to  be  widely  applicable  to  aeroelastic  problems  for  t urbomachinery 
cascades,  control  surface  freeplay,  airfoil/wing  flutter  and  LCO  (e.g.,  Refs  [11,  14,  15,  16]). 
As  evidenced  by  their  recent  studies,  under  AFOSR/STTR  contract  support,  the  Frequency- 
Domain  method  with  the  Harmonic  Balance  scheme  is  about  two-orders  of  magnitude  faster 
than  conventional  time-marching  methods. 

For  example,  the  3D  Time  Linearized  code  can  generate  the  needed  aerodynamic  model 
in  one  day  for  a  given  Mach  number  and  flutter  points  take  less  than  1  minute  per  Mach 
number  on  a  workstation  computer. 
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1.3.3  Frequency-Domain  POD/ROM  EigenMode  Approach 

Reduced-Order  Models  (ROM)  EigenMode 

The  eigenmodes  of  a  time-linearized  system,  which  may  be  thought  of  as  aerodynamic  states, 
were  computed  and  subsequently  used  to  construct  computationally  efficient,  reduced  order 
models  of  the  unsteady  flow  field.  An  important  advantage  of  the  eigenmode  based  re¬ 
duced  order  modeling  technique  is  that  once  the  eigenmode  information  has  been  computed, 
reduced  order  models  can  be  constructed  and  used  to  calculate  the  response  at  different 
frequencies  and  mode  shapes  for  almost  no  additional  computational  effort.  Furthermore, 
only  a  very  few  eigenmodes  or  states  (degrees  of  freedom)  need  to  be  retained  in  the  model 
to  accurately  predict  the  unsteady  aerodynamic  response,  making  the  method  ideally  suited 
for  rapid  flutter  calculation  and  active  control. 

The  POD  Technique 

The  Proper  Orthogonal  Decomposition  (POD)  technique,  also  known  as  Karhunen-Loeve 
(Ref.  [12])  expansions,  was  originally  introduced,  to  determine  and  model  coherent  structures 
in  turbulent  flow  fields.  Using  the  POD  approach,  one  examines  a  series  of  ”  snapshots”  of 
experimental  or  computational  data,  each  at  a  different  instant  in  time.  These  solution 
snapshots  are  used  to  form  a  small  and  more  compact  eigenvalue  problem  that  is  solved  to 
determine  a  set  of  optimal  basis  functions  for  representing  the  flow  field. 

Frequency-Domain  POD/ROM 

Hall,  Thomas,  and  Dowell  (Ref  [11])  developed  a  frequency-domain  form  of  the  POD  tech¬ 
nique,  and  applied  it  to  transonic  flows  about  airfoils.  In  particular,  they  used  a  time- 
linearized  CFD  analysis  to  compute  unsteady  small-disturbance  flow  solutions  for  vibrating 
airfoils  in  the  frequency  domain  over  a  range  of  frequencies.  Basis  vectors  were  then  ex¬ 
tracted  from  this  frequency-domain  data  set  using  the  POD  technique.  The  resulting  basis 
vectors  were  then  used  to  construct  low  degree  of  freedom  reduced-order  models  of  the  un¬ 
steady  flow.  Finally,  the  reduced-order  aerodynamic  model  was  combined  with  a  structural 
dynamic  model  resulting  in  a  compact,  but  accurate,  flutter  model. 

Harmonic-Balance  Method  for  Nonlinear  LCO  (Refs  [14,  15]) 

In  conjunction  with  the  frequency-domain  POD/ROM  method,  the  harmonic  balance  (HB) 
technique  can  be  used  in  the  treatment  of  nonlinear  LCO  problems  for  oscillating  wings  (with 
stores).  The  harmonic  balance  technique  is  used  to  recast  the  proposed  Euler  equations 
into  a  set  of  ’’higher-order”  equations  of  like  harmonics.  After  introducing  a  pseudo-time 
term  into  the  harmonic-balance  equations  so  that  they  may  be  solved  by  time  marching, 
these  equations  are  solved  by  conventional  CFD  techniques.  The  method  has  a  number 
of  advantages  over  the  more  conventional  time  domain  solutions.  Because  the  solutions  are 
computed  in  the  frequency  domain,  the  time-marching  algorithm  is  only  used  to  converge  the 
solution  to  steady  state.  Thus,  accelerating  techniques,  including  pseudo  time-marching  with 
multi-grid  acceleration  can  be  used.  Moreover,  for  problems  where  ’’engineering  accuracy” 
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is  required,  the  harmonic  balance  series  can  be  truncated  into  just  a  few  harmonics.  The 
result  is  that  the  proposed  harmonic  balance  method  is  potentially  two-orders  of  magnitude 
faster  than  conventional  time-marching  methods  for  determining  LCO. 

1.3.4  Merits  of  the  Frequency-Domain  Approach 

There  are  several  definite  advantages  of  using  the  frequency-domain  approach.  First,  the 
frequency-domain  harmonic  balance  method,  when  applied  to  high-level  equations,  can  be 
at  least  two-orders  of  magnitude  faster  than  the  nonlinear  time-domain  CFD  simulations. 
Also,  this  method  retains  essential  nonlinear  features  in  an  aeroelastic  system,  including  non¬ 
linear  structural  stiffness  and  damping  as  well  as  large  transonic  shock  excursions  including 
viscosity.  Second,  current  transonic  flutter  methods  using  time-domain/time-marching  CFD 
in  conjunction  with  various  versions  of  the  indicial  approach  (Refs  [17,  18])  are  very  tedious 
and  computationally  expensive  and  their  accuracy  usually  depends  on  the  indicial  motion 
imposed.  By  contrast,  the  frequency-domain  formulation  of  POD/ROM  directly  solves  for 
convenient  aerodynamic  mode  shapes  which  can  be  stored  and  repeatedly  applied  for  a 
large  range  of  frequencies;  hence,  its  a  much  more  efficient  and  accurate  method.  Third, 
the  frequency-domain  approach  with  an  eigenvalue  solution  method  is  a  familiar  practice 
to  the  structures/loads  and  flutter  engineers.  The  proposed  method  after  ZONA’s  further 
modification  is  expected  to  be  well  accepted  within  the  industrial  environment. 


1.4  The  ZONA/Duke  Team 

ZONA  Technology,  Inc.  (ZONA)  and  Duke  University  (Duke)  have  formed  a  strong  team 
with  a  comprehensive  background  to  handle  the  challenging  tasks  in  Phase  II.  Professor 
Hall  (P.I.)  and  Professor  Dowell,  one  of  the  world’s  leading  aeroelasticians,  at  Duke  are  the 
originators  of  the  proposed  ROM/POD  method,  which  has  gained  much  attention  in  the 
aerospace  community  in  recent  years. 

P.C.  Chen  (P.I.)  and  Danny  D.  Liu  at  ZONA  are  the  small  business  counterpart  of  this 
STTR/Phase  II.  ZONA  has  been  extensively  engaged  in  the  R&D  of  LCO  including  its 
prediction  methodology,  control  and  F-16  and  F-18  LCO  data  correlation.  ZONA  has  col¬ 
laborated  closely  with  Lockheed-Martin/LMTAS  and  the  Seek  Eagle  office  of  Eglin  AFB  on 
F-16  wing/store  LCO  investigations.  The  results  including  the  definition  of  a  NSD  model 
were  presented  in  two  AIAA/SDM  papers  (Refs  [7,  19]).  Under  a  recent  NAVAIR  contract, 
ZONA  also  worked  closely  with  team  member  Boeing/  St  Louis  on  a  successful  R&D  in  the 
reconfigurable  adaptive  control  of  the  F-18  LCO.  (Ref  [20]).  Supported  by  NASA/Langley, 
onging  development  of  the  CFD/CSD  interfacing  using  BEM  solver  (Refs  [21,  22])  will  pro¬ 
vide  a  new  spline  methodology  for  tightly  coupled  aerodynamic-structural  interaction  for  3D 
multiple  mode  LCO  studies  in  Phase  II. 

In  August  1999,  the  AFOSR  awarded  an  STTR  Phase  I  contract  to  the  ZONA  Technol¬ 
ogy  (ZONA)  team  to  develop  innovative  computational  aeroelasticity  methodologies  for  un- 
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derstanding,  predicting  and  controlling  critical  nonlinear  aerostructural  interaction  (e.g., 
LCO/flutter)  phenomena.  The  ZONA/Duke  team  has  successfully  accomplished  the  STTR 
Phase  II  contract. 


In  what  follows,  we  present  the  technical  content  of  the  Phase  II  final  report. 

Chapter  2  presents  a  comprehensive  LCO  and  flutter  study  of  a  nonlinear  inviscid  tran¬ 
sonic  unsteady  flow  over  a  NACA64A010A  airfoil  (Euler  Method).  Chapter  3  presents  a 
comprehensive  LCO  and  flutter  study  of  a  nonlinear  viscous  transonic  unsteady  flow  over 
a  NLR7301  supercritical  airfoil  (RANS  Method).  Both  studies  adopts  Duke  University’s 
innovative  Harmonic  Balance  method  in  the  frequency  domain.  As  a  counterpart  to  Chap¬ 
ter  3,  Chapter  4  presents  the  time-domain  simulation  of  the  nonlinear  viscous  flow  over 
the  same  NLR7301  airfoil  for  LCO  investigation  and  verification  (RANS  Method:  CFL3D 
code).  Chapter  5  presents  the  theoretical  background  of  Duke’s  Harmonic  Balance  method 
in  three-dimensions.  Chapter  6  summarizes  the  ZONA  effort  in  providing  the  structural 
modes  and  deforming  CFD  meshes  for  various  wing  planfrom  candidates  to  support  Duke’s 
LCO/flutter  investigation  using  Duke’s  Harmonic  Balance  in  three-dimensions.  Chapter  7 
presents  a  comprehensive  LCO  and  flutter  study  of  a  nonlinear  inviscid  transonic  unsteady 
flow  over  an  AGARD  445.6  wing  (3D  Euler  Method).  Chapter  8  presents  the  LCO  and  flut¬ 
ter  study  of  a  nonlinear  viscous  transonic  unsteady  flow  over  an  F-16  wing  including  effects 
due  to  stores  (3D  RANS  Method).  The  method  used  in  chapters  7  and  8  demonstrates  the 
significant  achievement  in  the  generalization  of  the  two-dimensional  Harmonic  Balance  (HB) 
into  a  three-dimensional  one  up  to  the  level  of  RANS  (Reynolds  Average  Navier-Stokes)  for 
its  LCO/flutter  application  to  a  complex  wing-store  system. 

This  Phase  II  final  report  intends  to  provide  a  comprehensive  description  of  our  approach 
and  results  obtained  in  achieving  the  Phase  II  technical  objectives  (see  1.2).  These  achieve¬ 
ments  provide  a  sound  foundation  upon  which  the  Phase  III  effort  for  technology  transfer 
for  a  production-ready  LCO  prediction  software  can  be  carried  out  thus  rendering  viable 
commercialization. 
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Chapter  2 

Transonic  Flutter  and  LCO  in 
Nonlinear  Inviscid  Two  Dimensional 
Flow:  Frequency  Domain 


Summary 

•  A  parametric  study  of  transonic  flutter  and  LCO  has  been  carried  out  for  a  NACA 
64A010A  airfoil  in  a  two-dimensional,  inviscid  (Euler)  flow.  Among  the  principal 
findings  are  that  the  LCO  may  be  either  benign  (stable)  with  the  LCO  existing  only 
beyond  the  flutter  boundary  or  dangerous  (unstable)  with  LCO  occurring  below  the 
flutter  boundary  if  a  sufficiently  large  disturbance  is  given  to  the  airfoil.  This  change 
from  benign  to  dangerous  LCO  occurs  when  the  Mach  number  is  varied,  for  example. 
Also  it  is  found  that  LCO  of  either  type  will  only  be  detectable  over  a  narrow  range  of 
Mach  number.  Outside  of  this  Mach  number  range  the  predicted  LCO  amplitudes  are 
so  large  as  to  be  indistinguishable  from  flutter. 

•  The  work  in  this  chapter  has  been  presented  at  the  AIAA/SDM  2002  Conference:  ”D.B. 
Kholodar,  J.P.  Thomas,  E.H.  Dowell,  and  K.C.  Hall,  A  Parametric  Study  of  Tran¬ 
sonic  Airfoil  Flutter  and  Limit  Cycle  Oscillation  Behavior,  AIAA  Paper  2002-1211, 
presented  at  the  AIAA/ASME/ASCE/AHS  SDM  Conference,  April  2002,  Denver,  CO; 
it  starts  with  the  following  abstract: 

•  Using  a  state  of  the  art  inviscid  computational  fluid  dynamic  harmonic  balance  aero¬ 
dynamic  Euler  based  code,  a  systematic,  parametric  investigation  is  presented  into  how 
the  several  structural  parameters  and  freestream  Mach  number  of  a  transonic  flow  affect 
the  flutter  characteristics  of  a  “typical”  two  degree- of -freedom  transonic  airfoil  config¬ 
uration.  The  computational  efficiency  of  the  time  linearized  option  of  the  harmonic 
balance  aerodynamic  model  allows  a  much  more  thorough  exploration  of  the  parameter 
range  than  has  been  possible  previously. 

2.1  Introduction 

Transonic  flow  flutter  and  limit  cycle  oscillations  (LCO)  are  of  significant  interest  in  wing 
and  aircraft  design.  The  large  expense  incurred  in  both  time-  and  frequency-domain  tran- 


8 


sonic  aerodynamic  computations  is  the  principal  obstacle  to  the  aeroelastician  in  obtaining 
a  deeper  understanding  of  these  phenomena  through  a  systematic  parameter  study. 

Reduced  order  modeling  (ROM)  techniques  have  been  developed  and  used  to  overcome  this 
obstacle  in  recent  work  on  this  subject.  For  a  general  review  of  the  latest  studies  involving 
ROM  based  methods  see  Ref.  [23].  Pade  approximants  of  transfer  functions  and  various 
other  rational  polynomial  curve  fitting  techniques  also  have  been  used  in  recent  years  for 
linear  flutter  analysis,  e.g.  see  Refs.  [24,  25,  26]. 

In  the  past  few  years  at  Duke  University  a  number  of  CFD  time  (dynamically)  linearized 
codes  have  been  developed  [27,  28]  and  converted  to  the  frequency  domain.  ROM  techniques 
have  then  been  applied  to  these  dynamically  linearized  CFD  codes  and  then  used  for  flutter 
analyses  that  are  very  computationally  efficient. 

Recently,  a  novel  nonlinear  harmonic  balance  (HB)  method  that  extends  the  frequency  do¬ 
main  CFD  models  to  the  fully  dynamically  nonlinear  range  has  been  developed.  This  method 
enables  one  to  model  efficiently  nonlinear  unsteady  aerodynamic  behavior  corresponding  to 
finite  amplitude  structural  motion  of  a  prescribed  frequency,  and  which  can  be  subsequently 
used  for  modeling  LCO  behavior  [27,  28,  29].  We  believe  these  two  methods,  ROM  based 
flutter  analysis  and  HB  based  LCO  modeling,  will  significantly  advance  the  aeroelastician’s 
capability  to  do  rapid  parametric  studies. 

In  the  current  study,  a  time  linerized  option  of  the  Euler  HB  model  is  used  to  capture 
the  effects  of  the  mean  position  of  the  shock  and  small  shock  motions  about  this  mean  po¬ 
sition  on  transonic  flutter.  The  shock  motion  is  assumed  linearly  proportional  to  the  airfoil 
motion  in  this  study. 

This  study  also  had  another  goal,  i.  e.  finding  a  (flutter)  boundary  of  neutrally  stable 
points  for  further  use  in  a  subsequent  Limit  Cycle  Oscillation  (LCO)  study  [29]  of  the  same 
model  with  large  shock  motion  aerodynamic  nonliiiearities  included. 

Results  for  the  flutter  boundary  using  the  present  methodology  for  the  famous  Isogai  typical 
section  case  have  been  discussed  in  Ref.  [30]  and  a  favorable  comparison  made  with  the 
results  of  other  investigators  who  used  potential  or  Euler  flow  models.  Further  comparisons 
with  results  of  others  from  theory  and  experiment  are  included  in  this  paper. 


2.2  Governing  Equations 

Consider  a  “typical”  two  degrees-of-freedom  (DOF)  airfoil  section  as  shown  in  Fig.  2.1.  The 
equations  of  motion  of  this  aeroelastic  system  can  be  written  in  the  form: 


mh  +  Saa  +  Khh  —  —L 
Sah  +  Iaa  +  KqPl  =  Mq.  a. 


(2.1) 
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Figure  2.1:  Sketch  of  Aeroelastic  Typical  Section. 


Here  the  left-hand  side  terms  represent  a  linear  structural  model  approximation  for  the 
plunge  and  pitch  coordinates.  The  right  hand  side  terms  represent  the  aerodynamic  load¬ 
ing  terms,  which  for  this  study  are  based  upon  the  harmonic  balance  approach  applied  to 
a  discrete  CFD  model  of  the  inviscid  Euler  equations.  Here  we  present  a  summary  of  the 
method.  For  a  more  detailed  description  see  Refs.  [27,  28]. 

In  integral  form,  the  unsteady  Euler  equations  may  be  written  as 

1 1 L vdD + L  (F(u)  -v%)dy-  L  (G(u)  -  uSk = °  (2-2) 

where  x  and  y  are  the  Cartesian  coordinates,  t  is  time,  and  D  is  a  deforming  control  volume 
bounded  by  the  control  surface  dD.  The  quantities  df/dt  and  dg/dt  are  the  x  and  y  com¬ 
ponents  of  the  velocity  of  the  control  surface  dD.  The  conservation  fluid  variables  U  and 
flux  vectors  (F,  G)  are  given  by 

U  =  {p,  pu ,  pv,  E}T 

F(U)  =  {pu,  pu2 +p,  puv,  (E +  p)u}T  (2.3) 

G(U)  =  {pv,  puv,  pv2  +p,  (E  +  p)v}T 

where  p  is  the  density,  u  and  v  are  the  velocity  components  in  the  x  and  y  directions  re¬ 
spectively,  p  is  the  pressure,  and  E  is  the  total  energy,  which  is  the  sum  of  the  internal  and 
kinetic  energy: 


Considering  a  calorically  perfect  gas,  the  system  of  equations  is  completed  via, 

p  =  (7  - 1)  {e  ~  |  K + v2] } 

2.3  Frequency  Domain  Equations 


(2.5) 


The  next  step  in  the  HB  development  is  to  consider  strictly  periodic  unsteady  flows  of  a 
fundamental  frequency  u.  The  flow  variables  and  flux  vectors  are  then  approximated  using 
truncated  Fourier  series  in  time  with  spatially  varying  coefficients.  For  instance, 


[[  Ud ^  U ne*™* 

hi)  n=-N 


and 


(2.6) 


i  (F(V)-V^-)dy-<f  (g(U)-u||')*:«  T  <2'7) 

JdD(t)  \  OXj  JdD(t)  \  oyj 

where  N  is  the  number  of  harmonics  used  in  the  Fourier  expansion.  The  time  derivative  of 
the  vector  of  conservation  fluid  variables  is  then 


M 

D(t) 


N 


U  dD  rj  inu  ^  U„e* 

n=—N 


inojt 


(2.8) 


Substituting  (2.7)  and  (2.8)  into  the  Euler  equations  (2.2),  multiplying  by  and  inte¬ 

grating  over  one  period  yields  a  system  of  equations  for  the  Fourier  coefficients 


AU  +  F  =  0 


(2.9) 


where  the  diagonal  matrix  A  and  vectors  U  and  F  are 


diagA  =  {-iN,  —  i(N  - 

U  =  {U-jv,  U(jv-i))  •••,  Un}t  (2.10) 

F  =  {F^r,F(jv-i),  ...,Fjv}t 


As  demonstrated  by  Hall  [27]  via  a  Fourier  transform  matrix,  E,  one  can  relate  the  Fourier 


coefficient  variables,  U,  to  the  solution  variables  at  (2 N  +  1)  discrete  time  levels  within  a 
period  of  motion.  Hall  also  showed  how  one  can  express  the  Fourier  Flux  Coefficients  in 
terms  of  time  domain  Flux  term  variables.  This  enables  the  HB  methodology  to  be  easily 
implemented  within  an  existing  steady  CFD  flow  solver  method.  If  one  defines 


and 


r  u(*o)  i 

//  U(to)dD  ) 

u  =  < 

Ditrj) 

i  i 

>  =  i 

...  ^ 

l  0(tK)  J 

//  UMd D 

(2.11) 


where 


then 


r  F(to)  \ 


>  = 


l  F^n)  ) 

[  §  [(F(U)  -  Ug)  \tody  -  (g(U)  -  Ug)  \todx] 


dD(t0) 


f  [(F(U)  -  Ug)  \t„dy  -  (g(U)  -  Ug)  | t„dx] 

^  dD(tN)  L  J 


I 


27rn 


(2  N  +  IV 


,  n  =  0, 1,  ...,2N 


U  =  EU,  F  =  EF 


(2.12) 

(2.13) 

(2.14) 


Substituting  Eqs.  (2.14)  into  the  equation  for  the  Fourier  coefficients,  Eq.  (2.9),  yields 

AEU  +  EF  =  0  (2.15) 

which  after  premultiplying  by  E-1  results  in 

DU  +  F  =  0  (2.16) 

where  D  =  E_1AE  is  the  spectral  source  term  operator  that  represents  and  approximates 
the  temporal  partial  derivative,  d/dt.  As  noted  in  Refs.  [27,  28],  to  solve  the  harmonic 
balance  equations,  one  can  add  a  pseudo  time  derivative  term  dtj/dr  to  Eq.  (2.16): 

+  DU  +  F  =  0  (2.17) 
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(a)  Overall  view 


(b)  Close-up 


Figure  2.2:  Depiction  of  the  97  x  97,  CFD  Grid  with  Outer  Boundary  Radius  of  Ten  Chord 
Lengths  around  a  NACA  64A010A  Airfoil. 


where  r  is  a  fictitious  time.  Since  only  a  “steady-state” ,  temporally  periodic  solution  is  de¬ 
sired,  local  time  stepping  and  multiple-grid  acceleration  techiniques  can  be  used  to  increase 
the  speed  of  computational  convergence.  The  above  described  procedures  allow  one  to  model 
nonlinear  unsteady  aerodynamic  of  a  prescribed  frequency  very  efficiently.  Again,  see  Refs. 
[27,  28]  for  additional  details. 


2.4  Computational  Model 

Fig.  2.2a  shows  the  “0”-type  computational  grid  used;  a  close-up  in  Fig.  2.2b  shows  the  grid 
cells  that  surround  a  symmetric  NACA  64A010A  airfoil.  The  mesh  consists  of  Nr  x  Nq  radial 
and  circumferential  nodes.  Flutter  results  presented  in  this  paper  are  computed  using  either 
65  x  65  or  97  x  97  grids.  However,  grids  of  25  x  25,  33  x  33  49  x  49  and  129  x  129  were  also 
considered  for  a  grid  convergence  study.  The  outer  boundary  radius  or  domain  radius,  R, 
was  chosen  to  be  either  5  or  10  chord  lengths  however  the  domains  of  R/c  =  15  and  R/c  =  20 
chord  lengths  were  also  considered.  The  dependence  of  the  calculated  aerodynamic  forces 
and  flutter  boundaries  on  the  grid  and  domain  size  was  investigated  to  ensure  convergence  of 
the  numerical  results.  A  slight  variant  of  the  standard  node-centered  Lax-Wendroff  scheme 
was  used  to  solve  the  Euler  equations.  See  Refs.  [31,  32]  for  further  details. 
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Figure  2.3:  Flutter  Speed  Index  (a)  and  Reduced  Frequency  (b)  vs  the  Squared  Value  of  the 
“Averaged”  Grid  Step  for  the  Meshes  N  x  N  of  N  =  129,  97,  65,  49,  33,  and  25. 
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0.50 


Figure  2.4:  Flutter  Speed  Index  (a)  and  Reduced  Frequency  (b)  vs  the  Domain  Radius. 
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2.5  Convergence  Study  Results 

The  dependence  of  the  calculated  aerodynamic  forces  on  the  grid  and  domain  size  was  investi¬ 
gated  for  Mach  numbers  of  M  =  0.7  and  M  =  0.8  for  various  values  of  the  reduced  frequency. 

In  the  case  of  the  grid  size  it  was  found  that  the  variation  of  aerodynamic  forces  is  monotonic 
and  “close  to  linear”  when  plotted  versus  a  squared  measure  of  grid  spacing.  This  is  what 
one  would  desire  since  the  CFD  code  used  is  a  second  order  model  in  the  spatial  variables. 
The  dependence  of  the  aerodynamic  forces  on  the  domain  size,  however,  is  not  monotonic. 
Shown  here  in  Figs.  2.3  and  2.4  are  the  dependence  of  the  flutter  boundaries  on  the  grid  and 
domain  size,  respectively.  In  Fig.  2.3  the  flutter  speed  index  and  reduced  frequency  appear 
to  be  nearly  linear  functions  of  the  grid  step.  Moreover,  noting  the  magnified  scales  of  the 
vertical  axes,  one  can  conclude  that  flutter  results  (for  these  Mach  numbers)  are  weakly  in¬ 
fluenced  by  the  grid  size.  The  flutter  results  in  Fig.  2.4  show  that  increasing  the  domain  size 
from  R/c  =  5  does  not  lead  to  substantial  changes  in  the  flutter  velocity  index  or  frequency. 
For  more  results  on  grid  and  domain  convergence,  see  Ref.  [33]. 
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Imaginary  Unsteady  Moment  Due  to  Plunge  Real  Unsteady  Moment  Due  to  Plunge 


Figure  2.5:  Real  and  Imaginary  Parts  of  the  Unsteady  Aerodynamic  Moment  due  to  Plunge 
Motion  vs  Reduced  Frequency  h  =  0.001. 


(b)  Ratio  of  Steady  Lift  to  Pitch  (c)  Ratio  of  Lift  Magnitude  to 

Angle,  ci /a,  Deduced  for  a  =  Pitch  Angle,  |cj/a|,  vs  a  for  dif- 

0.0001°  and  1°  vs  Mach  Number  ferent  reduced  frequencies  0  <  Q  < 

0.1  at  M  =  0.868. 


Figure  2.6:  Lift  Curve  Slope  Results. 
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2.6  Aerodynamic  Lift  and  Moment 

For  a  simple  harmonic  motion  of  the  airfoil, 

Y  =  heiuJt,  a  =  aeiw*. 

0 


(2.18) 


the  aerodynamic  forces  will  usually  be  periodic  in  time,  and  for  aeroelastic  analysis  the  first 
harmonic  will  be  dominant.  Thus  to  a  good  approximation, 


Me.  a.  = 


(2.19) 


In  general,  the  nondimensional  lift  and  moment  coefficients  are  expressed  as 

Ci  =  ci(u,  geometry,  M,h,  a) 

Cm  =  Cm(u,  geometry,  M,  h,  a) , 


(2.20) 


where  5/  and  Cm  are  nonlinear  functions  of  h  and  a.  However,  in  dynamically  or  time  lin¬ 
earized  aerodynamics  (used  in  flutter  analyses  per  se),  the  nondimensional  lift  and  moment 
coefficients  are  linearly  proportional  to  h  and  a: 


Ci  =  CiR (u,  geom,  M)h+cis(u,  geom,  M)a 
Cm=  geom,  M)h+ £„„(£,  geom,  M)a 


(2.21) 


For  reference  purposes,  the  case  of  zero  Mach  number  was  considered  as  well.  The  closed 
form  expressions  [33]  for  the  lift  and  moment  coefficients  due  to  the  pitch  and  plunge  motion 
were  obtained  from  the  classical  Theodorsen’s  formulas. 


Some  representative  results  for  linear  aerodynamic  behavior  computed  by  running  the  HB 
aerodynamic  solver  with  very  small  amplitude  motion  (h  =  0.001)  are  presented  in  Fig.  2.5, 
which  shows  the  real  and  imaginary  parts  of  the  aerodynamic  moment  due  to  the  plunge  mo¬ 
tion,  Cmh,  versus  reduced  frequency.  This  (symmetric  NACA  64A010A)  airfoil  has  a  strong 
shock  in  the  range  of  0.8  <  M  <  0.9.  For  subsonic  Mach  numbers  0  <  M  <  0.7,  Cmh  is 
relatively  constant  (on  the  reduced  frequency  interval  [0, 1]).  For  M  =  0.8  a  “wavy”  pattern 
is  very  prominent  in  Cmh.  This  “wavy”  pattern  reflects  the  slow  upstream  acoustic  waves 
typical  of  high  subsonic/transonic  flow  [34].  The  “wavy”  pattern  persists  in  the  aerodynamic 
moment  up  to  and  including  M  —  0.9  (the  Mach  numbers  of  0.82,  0.84,  0.86,  and  0.88  were 
also  studied,  but  the  results  axe  not  shown  in  Fig.  2.5).  At  higher  transonic  numbers,  from 
M  =  0.92  and  to  M  =  1.0,  the  shock  location  is  at  the  trailing  edge  of  the  chord,  the  flow 
over  the  airfoil  chord  is  supersonic,  and  the  “wavy”  pattern  is  no  longer  present.  (Not  plot¬ 
ted  here  are  curves  for  Mach  numbers  of  0.92,  0.94,  0.96,  and  0.98  as  they  correspond  very 
closely  to  M  =  1.0  results  shown  in  Fig.  2.5.)  The  oscillatory  behavior  of  Cmn  is  attributed  to 
the  tendency  of  the  fluid  dynamic  eigenspectrum  to  become  less  damped  for  transonic  Mach 
numbers.  See  Ref.  [23].  Previous  investigators  have  noted  that  computational  waviness  of 
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this  sort  may  be  due  to  inadequate  computational  grid  convergence,  Ref.  [35].  However,  it 
has  also  been  known  since  the  classical  results  of  Landahl  [36]  that  such  waviness  may  arise 
from  physical  effects.  The  grid  convergence  study  previously  discussed  confirmed  that  the 
waviness  in  the  present  study  is  physical  in  origin. 

Shown  in  Fig.  2.6a  is  a  plot  of  the  ratio  of  steady  lift  to  pitch  angle  versus  Mach  num¬ 
ber  (for  u  =  0)  as  obtained  from  this  Euler  flow  model.  For  reference  the  classical  results  as 
obtained  from  potential  thin  airfoil  theory,  dci/da  =  27t/a/1  —  M  at  a  =  0,  are  also  shown. 
The  effects  of  airfoil  profile  (thickness)  are  clearly  seen.  Fig.  2.6a  presents  the  Euler  flow 
model  result  for  C\  divided  by  a  (in  radians)  for  a  =  1°.  However,  it  is  found  that  for  the 
reduced  frequency  of  zero,  u  =  0,  Ci/a  is  an  erratic  function  of  Mach  number  in  the  range 
0.83  <  M  <  0.89  for  sufficiently  small  a  ,  see  Fig.  2.6b.  Note  this  is  true  even  though  for  a 
fixed  M,  ci/ a  is  a  smooth  function  of  a,  see  Fig.  2.6c.  However,  for  the  reduced  frequencies 
of  interest  for  the  flutter  boundaries  studied  here,  u>  >  0.1,  these  erratic  variations  in  cja 
are  not  present.  To  illustrate  this  the  case  of  Mach  number  M  =  0.868  is  considered.  As  one 
can  see  in  Fig.  2.6b,  at  this  value  (marked  by  the  vertical  line)  of  M  there  is  a  noticeable 
variation  in  the  values  of  ci/a  with  a  when  u;  =  0.  By  contrast,  as  shown  in  Fig.  2.6c  for 
M  =  0.868,  as  the  reduced  frequency  increases  above  Q  =  0.05,  the  ratio  of  lift  to  pitch  angle, 
Ci/ a,  is  virtually  constant  with  respect  to  a.  Similar  results  were  obtained  for  other  Mach 
numbers,  but  are  not  presented  here.  It  was  also  found  that  the  degree  of  numerical  smooth¬ 
ing  in  the  CFD  code  could  change  the  details  of  the  results  for  lift  at  very  small  angles  and 
very  low  reduced  frequencies  in  the  Mach  number  range  where  a  rapid  variation  in  lift  occurs. 

Finally,  a  comparison  of  aerodynamic  transfer  functions  for  this  (NACA  64A010A)  air¬ 
foil  has  been  performed  for  M  =  0.8.  See  Fig.  2.7.  Lift  qa  and  moment  Cma  coefficients 
(both  per  radian)  obtained  by  the  present  Euler  code  are  compared  to  experimental  values 
of  Davis  and  Malcom[37],  other  Euler  calculations  performed  by  Magnus[38]  and  Bendiksen 
and  Kousen[39].  Also  shown  are  the  results  of  IJeda  and  Dowell[40]  obtained  by  a  describing 
function  method  based  on  LTRAN2  solver  for  transonic  small-disturbance  potential  theory 
as  well  as  similar  results  by  Isogai[41].  The  results  of  Refs.  [37,  38,  39,  40,  41]  shown  in 
Fig.  2.7  were  taken  from  Figs.  3  and  4  of  Ref.  [39].  The  results  obtained  by  the  present 
Euler  code  were  computed  using  one  harmonic  in  the  aerodynamic  solution.  No  significant 
variations  in  lift  and  moment  coefficients  for  the  fundamental  harmonic  were  observed  when 
calculations  were  performed  using  two  harmonics. 

The  several  theoretical  results  shown  in  Fig.  2.7  for  the  lift  and  moment  as  a  function 
of  reduced  frequency  are  on  the  whole  consistent.  However  there  is  some  significant  varia¬ 
tion  among  them  especially  for  phase  angle,  but  less  variation  for  magnitude.  The  trends 
predicted  by  the  present  theory  and  those  from  experiment  are  in  encouragingly  good  agree¬ 
ment.  Previous  theoretical  work  provided  data  at  relatively  few  reduced  frequencies.  The 
computational  efficiency  of  the  present  theoretical  method  has  allowed  a  more  complete  data 
set  to  be  obtained. 
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(a)  Unsteady  Lift  due  to  Pitch  Magnitude,  \ci-\ 


(b)  Unsteady  Lift  due  to  Pitch  Phase  Angle 


(c)  Unsteady  Moment  due  to  Pitch  Magnitude,  (d)  Unsteady  Moment  due  to  Pitch  Phase  Angle 

I 


Figure  2.7:  Lift  and  Leading-Edge  Moment  Coefficients  (both  per  Radian)  as  Functions  of 
Reduced  Frequency  k  =  yr-  due  to  Pitching  ±1°  at  the  quarter-chord  for  M  =  0.8. 

Uoo 
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Mach  number,  M 

Mach  number,  M 

(a)  Flutter  Speed  Index, 

(b)  Flutter  Frequency  Ratio, 
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Mach  number,  M 

(c)  Flutter  Reduced  Frequency,^/ 
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Mach  number,  M 

(d)  Amplitude  of  Flutter  (Aeroelastic)  Eigen 
vector,  \hf\/ar 


Figure  2.8:  Flutter  Behavior  as  a  Function  of  the  Mass  Ratio  for  ^  =  0.5. 
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2.6.1  Mass  Ratio  Effect 


Dependence  of  the  flutter  boundaries  (flutter  speed  index,  flutter  frequency  and  eigenvector) 
on  the  mass  ratio,  p,  was  investigated  first.  Two  cases  for  the  ratio  of  uncoupled  frequencies 
are  presented  here:  uJh/ua  =  0.5  (Fig.  2.8)  and  uJh/u<x  =  0.8  (Fig.  2.9). 

^  =  0.5: 

For  uih/uia  =  0.5,  the  flutter  speed  index  (Fig.  2.8a)  is  seen  to  be  weakly  dependent  on  p 
for  subsonic  Mach  numbers.  For  transonic  Mach  numbers  of  M  =  0.8  and  M  =  0.9  the 
flutter  speed  index  decreases  with  p  at  a  much  higher  rate.  For  M  =  1.0  the  flutter  speed 
index  becomes  very  large  (well  beyond  the  scale  of  the  graph).  Note  that  the  flutter  reduced 
frequency  (Fig.  2.8c)  for  M  =  0.9  is  distinctly  higher  than  for  other  Mach  numbers.  The 
flutter  frequency  (Fig.  2.8b)  and  flutter  eigenvector  (Fig.  2.8d)  are  also  less  sensitive  to  mass 
ratio  at  low  subsonic  Mach  numbers  with  a  more  significant  dependency  occurring  at  higher 
subsonic  and  transonic  Mach  numbers.  Coupled  in  vacuuo  natural  frequencies,  U\  and 
were  also  computed  for  easy  reference  and  as  can  be  seen  in  Fig.  2.8b,  the  flutter  motion 
is  pitch  dominated  for  M  =  0.9.  One  may  also  see  this  from  Fig.  2.8d,  where  coupled  in 
vacuuo  eigenvectors,  {h/a) i  and  {h/a) 2,  are  indicated.  The  flutter  eigenvector  for  M  =  0.9 
is  essentially  the  same  as  {h/a)  1,  which  corresponds  to  the  pitch  dominated  mode.  This  is 
an  example  of  single  DOF  flutter,  but  note  the  critical  aeroelastic  mode  is  a  coupled  natural 
mode,  albeit  one  that  is  pitch  dominated. 


For  0Jh/uJa  =  0.8,  the  flutter  speed  index  (Fig.  2.9a)  is  even  more  weakly  dependent  on  p  for 
subsonic  Mach  numbers.  In  Fig.  2.9c,  the  flutter  speed  index  behavior  for  transonic  Mach 
numbers  of  M  =  0.8,  M  =  0.9  and  M  =  1.0  is  shown.  As  in  the  case  of  u)h/ua  =  0.5,  the 
flutter  speed  index  for  the  transonic  Mach  numbers  is  more  sensitive  to  p  than  for  the  lower, 
subsonic  Mach  numbers.  Moreover,  for  this  frequency  ratio,  when  M  =  0.9,  multiple  flutter 
velocities  occur  for  mass  ratios  p  >  80.  To  determine  the  stability  of  the  regions  created  by 
such  multiple  branches  of  flutter  boundaries  a  root-loci  analysis  for  fixed  values  of  Uh/uja 
and  p  can  be  employed  [42,  ?,  43]  using  reduced  order  aerodynamic  models.  However  in 
this  study  only  flutter  boundaries  were  determined.  The  flutter  frequency  (Fig.  2.9b)  and 
flutter  eigenvector  (Fig.  2.9d)  are  distributed  over  the  full  range  between  coupled  natural 
frequencies  and  eigenvectors  and  none  of  Mach  number  cases  appears  to  be  primarily  pitch 
or  plunge  dominated. 


2.6.2  Ratio  of  Uncoupled  Natural  Frequencies  Effect 

Dependence  of  the  flutter  boundary  on  the  ratio  of  uncoupled  natural  frequencies, 
is  next  investigated,  and  two  mass  ratios  are  considered:  p  =  25  (Fig.  2.10  and  2.11)  and 
p  =  100  (Fig.  2.12). 
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vf/Vfi 


(a)  Flutter  Speed  Index, 


(c)  Flutter  Reduced  Frequency,  u ;/  , 


Mach  number,  M 


(b)  Flutter  Frequency  Ratio,  ^ 


0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0 

Mach  number,  M 

(d)  Amplitude  of  Flutter  (Aeroelastic)  Eigen¬ 
vector,  \hf\/ar 


Figure  2.10:  Flutter  Behavior  as  a  Function  of  the  Ratio  of  Uncoupled  Frequencies  for 
/i  —  25. 
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Qiop*i 


Figure  2.11:  Flutter  Reduced  Frequency,  Uf,  versus  a)  £■§■  and  b)  the  coefficient  a2  of  a2(^ 
,.,2  “ 

+  a0  =  Ofor  =  25. 
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vf/Vji  vf/Vn 


(a)  Flutter  Speed  Index, 


(b)  Flutter  Frequency  Ratio, 


(c)  Flutter  Speed  Index,  ^ 


(d)  Amplitude  of  Flutter  (Aeroelastic)  Eigen¬ 
vector,  \hf\/ar 


Figure  2.12:  Flutter  Behavior  as  a  Function  of  the  Ratio  of  Uncoupled  Frequencies  for 
fi  =  100. 
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//  =  25: 


For  n  =  25,  the  flutter  speed  index  (Fig.  2.10a)  has  a  minimum  near  ojh/ua  =  1.  For  typical 
airfoils,  oJh/u>a  is  usually  between  0.2  and  0.8.  In  that  range,  the  flutter  speed  index  decreases 
steadily  and  is  almost  constant  for  Mach  numbers  less  than  0.7.  However,  as  can  be  seen 
and  is  well  known,  the  flutter  speed  index  becomes  very  sensitive  to  the  Mach  number  for 
high  subsonic  or  transonic  conditions. 

For  a  better  understanding  of  how  the  flutter  speed  index  results  are  obtained  one  can 
show  that  the  flutter  reduced  velocity  may  be  expressed  as 


(2.22) 


where  di,  d2,  and  are  coefficients  dependent  on  the  reduced  frequency.  Eq.  (2.22)  is 
obtained  from  the  imaginary  part  of  Eq.  (??).  Substitution  of  Eq.  (2.22)  into  the  real  part 
of  Eq.  (??),  leads  to  a  quadratic  equation  in  “l/ul, 


with  the  roots  given  by 


(2.23) 


_  — Qi  ±  \/af  -  4a2Qo 
\u>%J  1,2  2  a2  ’ 


where  the  coefficients  a»  may  also  change  the  sign  for  different  O.  Obviously,  only  positive 
real  values  of  in  Eq.  (2.24)  have  physical  meaning.  Moreover,  only  such  values  of 

Wfc/w*  are  kept  for  which  Vf  in  Eq.  (2.22)  is  positive.  Real  roots  of  Eq.  (2.24)  are  shown  in 
Fig.  2.11a.  When  the  coefficient  a2  passes  through  zero,  branches  of  (*>%/<*>%  assymptotically 
go  to  ±oo.  The  dependence  of  a2  with  respect  to  Of  is  shown  in  Fig.  2.11b.  For  example, 
in  Fig.  2.10c  for  M  =  0.9,  a2  goes  through  zero  at  Of  ~  0.7,  which  is  when  the  branches  of 
^h/^a  S°  to  ioo  in  Fig.  2.11a. 

Note  that  the  difference  between  the  two  coupled  natural  frequencies  (u^  and  u2  as  in¬ 
dicated  in  Fig.  2.10b)  is  the  smallest  near  ujh/ua  =  1.  This  is  likely  responsible  for  the 
minimum  in  the  flutter  speed  index  in  the  vicinity  of  that  point:  the  closer  the  coupled 
natural  frequencies  are  initially,  the  more  readily,  when  increasing  velocity,  “coalescence”  or 
“merging  frequency”  flutter  can  occur  (see  p.  112  of  Ref.  [34]). 


//  =  100: 

For  /x  =  100,  the  flutter  speed  index  (Fig.  2.12a  and  c)  again  has  the  tendency  to  have  a 
minimum  near  uih/oja  =  1,  and  especially  so  for  subsonic  Mach  numbers  (Fig.  2.12a).  As  in 
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the  case  for  n  =  25,  the  flutter  speed  index  is  nearly  constant  with  M  for  subsonic  condi¬ 
tions.  At  higher  transonic  Mach  numbers  (Fig.  2.12c),  the  flutter  speed  index  becomes  very 
sensitive  to  Mach  number. 

This  behavior  is  associated  with  the  “transonic  flutter  dip”.  Also  note  that  the  flutter 
speed  index  can  be  a  multivalued  function  of  u)h/ua  for  M  =  0.8  and  M  =  0.9. 

Note  that  (Fig.  2.12b)  the  flutter  frequency  is  very  close  to  the  coupled  natural  pitch  fre¬ 
quency,  u\/uja,  for  M  =  0.9  (pitch  dominated  motion).  This  is  indicative  of  “single-degree- 
of-freedom”  flutter  (see  p.  113  of  Ref.  [34]). 

2.6.3  Mach  Number  Effect 

The  dependence  of  the  flutter  boundary  on  Mach  number  (see  Figs.  2.13  and  2.14)  for  the 
cases  of  fj,  =  25,  100  and  ix>h/u>a  =  0.5,  0.8  is  presented  next. 

^  =  0.5: 

Fig.  2.13  shows  the  flutter  boundary  behavior  versus  Mach  number  for  u>h/uJa  =  0.5.  Squares 
represent  the  case  of  //  =  25  and  circles  the  case  of  n  =  100.  Note  a  rather  weak  dependence 
of  the  flutter  boundaries  on  the  mass  ratio. 

As  the  Mach  number  is  increased  past  M  =  0.8,  a  sharp  drop  in  the  flutter  reduced  frequency 
occurs.  That  is  accompanied  by  a  sharp  increase  in  the  flutter  speed  index.  As  the  Mach 
number  is  increased  further,  another  branch  of  flutter  points  appears.  These  low  flutter 
speed  index  values  correspond  to  the  “transonic  flutter  dip” . 

For  an  explanation  to  the  abrupt  change  in  flutter  mode  at  M  =  0.82  and  M  =  0.91 
one  may  ask  whether  the  steady  flow  shock  position  is  changing  and  where  the  shock  loca¬ 
tion  is  relative  to  the  elastic  axis  and  mode  nodes. 

The  steady  flow  pressure  contours  for  M  =  0.81  and  M  =  0.82  did  not  show  any  significant 
changes  in  the  pressure  distribution.  However  the  situation  is  different  for  Mach  numbers 
above  M  =  0.9.  Here  the  location  of  the  shock  moves  to  the  trailing  edge  as  the  Mach 
number  is  increased  above  M  =  0.9.  For  uJh/uja  =  0.8,  the  flutter  behavior  also  changes  in 
the  viscinity  of  these  Mach  numbers.  Note  that  in  the  Mach  number  range  0.82  <  M  <  0.92 
(Fig.2.15)  the  position  of  the  shock  on  the  airfoil  moves  appreciably.  To  develop  an  im¬ 
proved  understanding  of  the  behavior  of  the  flutter  boundaries  one  could  perform  a  root 
locus  analysis  which  would  enable  one  to  view  the  migration  of  the  aeroelastic  eigenvalues 
as  a  function  of  the  nondimensional  airspeed  for  each  Mach  number.  Representative  root 
locus  analyses  for  transonic  Mach  numbers  of  similar  airfoil  models  can  be  found  in  Refs. 
[42,  30,  43].  Root  loci  however  have  not  been  part  of  the  present  study. 

The  reduced  flutter  frequency  behavior  is  shown  in  Fig.  2.13c.  For  subsonic  and  low  transonic 
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vf/Vn 


(a)  Flutter  Speed  Index, 


(c)  Flutter  Reduced  Frequency,  a?/ 


(b)  Flutter  Frequency  Ratio,  (^L 


(d)  Amplitude  of  Flutter  (Aeroelastic)  Eigen¬ 
vector,  \hf\/ar 


Figure  2.13:  Flutter  Behavior  as  a  Function  of  Mach  Number  for  ^  =  0.5  and  a  —  25,  100. 
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Vf/V|A 


(a)  Flutter  Speed  Index, 


(c)  Flutter  Reduced  Frequency,  tof 


(b)  Flutter  Frequency  Ratio, 


(d)  Amplitude  of  Flutter  (Aeroelastic)  Eigen¬ 
vector,  \hf\/ar 


Figure  2.14:  Flutter  Behavior  as  a  Function  of  Mach  number  for  ^  =  0.8  and  a  =  25,  100. 
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(a)  M  =  0.82  (b)  M  =  0.92 


Figure  2.15:  Steady  Flow  Pressure  Contours. 

Mach  numbers  (0  <  M  <  0.8),  Of  is  much  higher  for  the  lower  mass  ratio  (Of  «  0.4)  than 
for  the  higher  mass  ratio  (Of  «  0.2).  At  high  transonic  Mach  numbers  (0.82  <  M  <  0.92), 
the  flutter  reduced  frequency  is  about  the  same  for  either  mass  ratio. 

The  flutter  (aeroelastic)  eigenvector  is  shown  in  Fig.  2.13d  and  provides  a  possible  explana¬ 
tion  for  the  low  flutter  speed  index  values  for  the  transonic  dip  Mach  number  region.  In  this 
range,  the  flutter  mode  becomes  pitch  dominated  for  both  values  of  mass  ratio.  The  flutter 
mode  is  markedly  different  from  that  for  M  <  0.82,  for  example. 

^  =  0.8: 

The  flutter  speed  index  for  Uh/ua  =  0.8  is  shown  in  Fig.  2.14a.  For  this  rather  high  value 
of  pitch  to  plunge  frequency  ratio,  there  is  no  transonic  dip  in  the  flutter  speed  index  at 
the  lower  mass  ratio.  Moreover  the  system  appears  to  be  more  stable  at  the  high  transonic 
Mach  numbers.  For  the  higher  mass  ratio  as  the  Mach  number  increases,  the  flutter  velocity 
index  reaches  its  minimum  in  the  range  0.78  <  M  <  0.88.  The  flutter  mode  then  appears 
to  change  twice  for  even  higher  Mach  numbers.  Contrary  to  the  results  for  Uh/ua  =  0.5, 
here  the  flutter  frequency  shown  in  Fig.  2.14b  is  not  very  close  to  either  the  pitch  or  plunge 
coupled  natural  frequency.  So  the  flutter  mode  (Fig.  2.14d)  is  not  a  single  degree  of  freedom 
flutter  motion,  but  rather  a  combined  pitch-plunge  motion. 
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Conclusions 

Using  a  state  of  the  art  Euler  based  time  linearized  aerodynamic  code,  an  investigation  is 
presented  of  how  structural  and  aerodynamic  parameters  including  freestream  (transonic) 
Mach  number  affect  flutter  and  LCO  characteristics  of  a  “typical”  two  degree-of-freedom 
airfoil  configuration.  Convergence  of  the  CFD  code  with  respect  to  grid  spacing  and  com¬ 
putational  domain  size  was  also  considered. 

The  following  conclusions  have  been  drawn: 

I  When  examining  the  dependence  of  the  flutter  on  the  mass  ratio,  /x,  it  was  found  that 
the  flutter  variables  (i.e.  reduced  velocity,  frequency  and  structural  eigenmode)  vary 
little  with  /x,  especially  for  subsonic  Mach  numbers. 

II  As  expected,  a  study  of  the  effect  of  the  ratio  of  uncoupled  natural  frequencies,  w/,/wQ, 
determined  that  flutter  reduced  velocities  have  a  minimum  near  u>h/oJa  «  1.  Multiple 
values  of  the  flutter  velocity  can  occur  at  high  transonic  Mach  numbers  for  some  u>h/ ua. 

III  Most  significantly  perhaps,  it  was  demonstrated  that  flutter  solutions  are  very  sensitive 
to  Mach  number  especially  in  the  transonic  range.  Depending  on  the  frequency  and 
the  mass  ratios,  there  may  or  may  not  be  sudden  and  significant  changes  (e.g.  the 
“transonic  dip”)  in  the  flutter  reduced  velocity  and  the  type  of  flutter  motion  as  the 
Mach  number  is  varied. 

IV  Finally  it  is  noted  that  viscous  effects  may  be  important.lt  has  been  our  experience 
that  when  the  aerodynamic  flow  is  attached  an  inviscid  analysis  is  adequate,  but  for 
separated  flows  a  viscous  analysis  is  required. 

An  extensive  parameter  analyses  of  an  airfoil  aeroelastic  configuration  has  been  achieved 
using  a  highly  efficient  time  linearized  CFD  computational  technique.  Correlation  with 
experiment  remains  an  open  challenge,  however  a  comparison  of  results  from  our  CFD  model 
with  results  from  an  unsteady  aerodynamic  experiment  by  Davis  and  Malcom[37]  for  lift  and 
moment  shows  encouraging  agreement. 


Nomenclature 


a 

b,c 

Cl, 

Clfr  ,  Qa 


CmfliCm& 

e 

h,a 


nondimensional  location  of  airfoil  elastic  axis,  a  =  e/b 
semi-chord  and  chord,  respectively 

first  harmonic  coefficient  of  lift  and  moment  about  elastic  axis,  respectively 
first  harmonic  coefficient  of  lift  due  to  plunge  and  pitch  motions,  respectively 
Note:  Qa  =  ci/a  for  h  =  0,  as  a  —*  0  cla  =  dci/da\a  =  0. 

Similar  definitions  hold  for  cin,Cmh  and  Cm-. 

first  harmonic  coefficient  of  moment  due  to  plunge  and  pitch  motions,  respectively 
location  of  airfoil  elastic  axis,  measured  positive  aft  of  airfoil  midchord 
airfoil  plunge  and  pitch  degree-of-freedom,  respectively 
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h,  a  = 

havr  ~ 

la 

L,Me.  a. = 
Kh,Ka  = 
k 
M 

m  = 

H  = 

N 

Nr>N&  = 
R 

ra  = 

Sa 

Uoo  = 

V 

v,uu  = 

L 0,0)  = 

~~ 

CJi,  a>2  = 

Xa  = 


first  harmonic  amplitude  of  plunge  nondimesionalized  by  semi-chord,  h  =  h/b, 
and  pitch  motion,  respectively 

“averaged”  grid  step  in  the  radial  direction,  havr  =  (R/c  —  1/2 )/(Nr  —  1) 
second  moment  of  inertia  about  elastic  axis 
lift  and  moment  about  elastic  axis,  respectively 

airfoil  plunge  stiffness  and  torsional  stiffness  about  elastic  axis,  respectively 

reduced  frequency  based  on  airfoil  semi-chord,  k  =  ub/U0 0 

freestream  Mach  number 

airfoil  sectional  mass 

mass  ratio,  n  =  m/'Kp00b2 

number  of  harmonics  used  in  harmonic  balance  solver 

number  of  grid  points  in  radial  and  circumferential  direction,  respectively 

radius  of  computational  domain 

radius  of  gyration  of  airfoil  about  elastic  axis,  =  Ia/mb 2 

first  moment  of  inertia  about  elastic  axis 

freestream  velocity 

reduced  velocity,  V  =  Uoo/uab 

flutter  speed  index,  Vf/y/Ji  =  Utx,f/y/jiuab 

frequency  and  reduced  frequency  based  on  airfoil  chord,  u  =  uc/Uqo 
uncoupled  natural  frequencies  of  pitch  and  plunge  DOF 
coupled  structural  natural  frequencies 
airfoil  static  unbalance,  xa  =  Sa/mb 


Subscripts 

/  =  flutter  onset  condition 
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Chapter  3 

Transonic  Flutter  and  LCO  in 
Nonlinear  Viscous  Two  Dimensional 
Flow:  Frequency  Domain 


Summary 

•  As  an  improvement  to  the  inviscid  model  presented  in  the  previous  chapter,  the  ef¬ 
fects  of  viscosity  have  been  added  and  initial  results  obtained  for  LCO.  Comparisons 
have  been  made  with  the  results  of  more  computationally  intensive  time  marching  CFD 
solutions  and  also  with  the  results  of  experiments  done  at  DLR  in  Gottingen  on  a  super¬ 
critical  airfoil.  The  HB  and  time  marching  results  are  in  good  agreement,  but  because 
of  the  greater  computational  efficiency  of  the  HB  method  a  much  broader  range  of  pa¬ 
rameter  values  has  been  considered  which  allows  a  more  complete  comparison  with  the 
experimental  data.  Based  upon  the  calculations  done  to  date,  it  appears  that  the  effects 
of  viscosity  on  LCO  are  most  pronounced  when  flow  separation  occurs  as  is  the  case 
for  the  DLR  study  of  a  supercritical  airfoil.  On  the  other  hand,  for  the  attached  flows 
that  occur  for  a  conventional  airfoil  (NACA  64 AO  10)  at  small  mean  angle  of  attack, 
the  results  from  the  inviscid  and  viscous  flow  calculations  are  more  closely  correspond. 

•  The  work  in  this  chapter  has  been  presented  at  the  AIAA/SDM  2002  Conference:  ”J.P. 
Thomas,  E.H.  Dowell,  and  K.C.  Hall,  Modeling  of  Viscous  Transonic  Limit  Cycle  Os¬ 
cillation  Behavior  Using  a  Harmonic  Balance  Approach,  AIAA  Paper  2002-1414,  pre¬ 
sented  at  the  AIAA/ASEM/ASCE/AHS  SDM  Conference,  April  2002,  Denver,  CO.”, 
it  starts  with  the  following  abstract: 


•  Presented  is  a  harmonic-balance  computational  fluid  dynamic  approach  for  modeling 
limit  cycle  oscillation  behavior  of  aeroelastic  airfoil  configurations  in  a  viscous  transonic 
flow.  For  the  NLR  7301  airfoil  configuration  studied,  accounting  for  viscous  effects  is 
shown  to  significantly  influence  computed  limit  cycle  oscillation  trends  when  compared 
to  an  inviscid  analysis.  Methodology  for  accounting  for  changes  in  mean  angle- of- attack 
during  limit  cycle  oscillation  is  also  developed. 
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3.1  Introduction 


A  novel  approach  for  modeling  limit  cycle  oscillation  (LCO)  behavior  of  aeroelastic  airfoil 
configurations  in  viscous  transonic  flows  is  presented.  The  method  is  based  on  a  harmonic 
balance  (HB)  flow  solver  technique  for  efficiently  modeling  nonlinear  unsteady  aerodynamics 
(see  Hall  et.  al.  [44]  and  Thomas  et.  al.  [45]).  The  primary  objective  of  the  current  study 
is  to  assess  the  capability  of  the  HB/LCO  solution  methodology  to  model  transonic  viscous 
flow  phenomena  such  as  shock  induced  boundary  layer  separation,  and  to  determine  if  such 
effects  have  an  influence  on  LCO  behavior.  Recent  two-dimensional  experimental  LCO  in¬ 
vestigations  by  Schewe  et.  al.  [46,  47]  and  Knipfer  and  Schewe  [48]  have  also  motivated  the 
present  research.  In  addition  to  viscous  effects,  and  unlike  the  NACA  64A010  airfoil  config¬ 
uration  studied  in  Thomas  et.  al.  [45],  we  now  also  consider  nonsymmetric  airfoil  sections 
at  nonzero  angles-of-attack.  As  will  be  shown  in  the  following,  this  requires  the  solution  of 
one  additional  equation  for  the  mean  (i.e.  zeroth  harmonic)  angle-of-attack  of  the  LCO. 

The  primary  motivation  for  developing  the  HB/LCO  solution  methodology  has  been  to 
construct  an  efficient  computational  procedure  for  modeling  LCO  behavior.  Figure  3.1a 
illustrates  a  ”  typical”  flutter  onset  boundary  in  the  reduced  velocity  versus  Mach  number 
plane,  including  the  commonly  observed  ’’flutter  speed  dip”.  As  one  moves  from  point  A, 
a  reduced  velocity  below  the  flutter  onset  condition,  to  point  B,  a  reduced  velocity  above 
the  flutter  onset  condition,  LCO  may  be  encountered.  Figure  3.1b  illustrates  three  possible 
types  of  LCO  behavior  one  might  observe. 

Curve  I  represents  a  stable  LCO  behavior  with  no  LCO  occurring  below  the  linear  flutter 
speed.  This  type  of  LCO  behavior  is  also  sometimes  referred  to  as  ’’soft  flutter”  or  ’’benign 
LCO”  per  se  in  that  the  nonlinear  effects  help  to  reduce  the  amplitude  of  the  oscillations.  Of 
course,  structural  failure  may  eventually  be  an  issue  if  the  amplitude  of  the  LCO  becomes 
too  large.  For  the  type  of  LCO  behavior  illustrated  by  Curve  II,  as  soon  as  one  reaches  the 
flutter  onset  condition,  any  further  increase  in  reduced  velocity  immediately  results  in  very 
large  LCO  motion.  Finally,  curve  III  illustrates  unstable  LCO  behavior.  A  configuration 
exhibiting  this  very  dangerous  type  of  LCO  trend  is  suseptible  to  ’’explosive”  LCO  whereby 
a  disturbance  of  sufficient  magnitude  may  be  capable  of  triggering  LCO  at  reduced  velocities 
even  below  the  flutter  onset  condition.  As  will  be  shown,  the  HB/LCO  solution  methodology 
is  a  computationally  efficient  technique  for  determining  such  LCO  trends. 


3.2  Theoretical  Development 

3.2.1  Fluid  Dynamic  Model 

For  the  model  configuration  studied  in  this  investigation,  we  consider  both  viscous  and  invis- 
cid  fluid  dynamic  models.  The  harmonic  balance  flow  solver  technique  provides  an  efficient 
method  for  modeling  nonlinear  unsteady  aerodynamic  effects  for  finite  amplitude  motions  of 
a  prescribed  frequency. 
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As  is  discussed  in  Hall  et.  al.  [44]  and  Thomas  et.  al.  [45],  the  HB  method  is  imple¬ 
mented  within  the  framework  of  a  conventional  CFD  solver.  The  vector  of  unsteady  flow 
conservation  variables  U(xi,t)  at  each  computational  mesh  point  x*  is  approximated  in  a 
truncated  Fourier  expansion  as 


nh 

U (3.1) 

71= — Njj 

where  Un(xj)  is  the  harmonic  coefficient,  and  NH  is  number  of  harmonics  used  in  the 
expansion.  McMullen  et.  al.  [49,  50]  are  also  investigating  a  similar  expansion  technique  for 
unsteady  Euler  and  Navier-Stokes  flows. 

Normally,  the  harmonic  balance  development  proceeds  by  first  substituting  Eq.  3.1  into 
the  governing  flow  equations.  One  then  subsequently  goes  through  the  process  of  ”  balanc¬ 
ing”  all  the  resulting  terms  proportional  to  ejnujt  for  each  n  (—NH  <n<  NH).  This  in  turn 
yields  2NH  +  1  equations  for  the  2 NH  +  1  harmonic  coefficients  Ura.  This  straight-forward 
approach  to  the  harmonic  balance  formulation  is  however  typically  difficult  to  implement  for 
complex  systems  of  equations  such  as  those  arising  from  Euler  and  Navier-Stokes  flows. 

Hall  et.  al.  [44]  however  recently  devised  an  alternative  approach  to  the  harmonic  balance 
derivation  whereby  the  method  is  formulated  in  terms  of  time-domain  variables.  That  is, 
instead  of  working  in  terms  of  the  Fourier  coefficient  variables  Un(xj),  one  instead  considers 
as  dependent  variables,  flow  solutions  stored  at  2NH  +  1  equally  spaced  ’’sub-time”  levels 
(U(xi,  tn))  over  the  period  of  one  cycle  of  motion.  The  Fourier  and  time-domain  variables 
are  in  fact  related  to  one  another  via  a  constant  Fourier  transformation  matrix. 

Working  in  terms  of  ”  sub- time”  level  variables  circumvents  the  necessity  of  having  to  go 
though  the  ’’balancing”  step  of  the  HB  method,  and  in  fact,  it  allows  one  to  easily  formulate 
the  harmonic  balance  technique  within  the  framework  of  an  existing  steady  CFD  solver  (See 
Hall  et.  al.  [44]  and  Thomas  et.  al.  [45]  for  further  details).  For  the  results  presented  here, 
the  CFD  method  is  a  variant  of  standard  Lax-Wendroff  scheme  [51,  52]  in  conjunction  with 
the  one-equation  turbulence  model  of  Spalart  and  Allmaras  [53] . 

3.2.2  Two-Dimensional  Airfoil  Aeroelastic  Model 

The  frequency  domain  form  of  the  unsteady  aeroelastic  equations  governing  the  ’’typical” 
pitch/plunge  airfoil  configuration  (Fig.  3.2)  can  be  written  as 

V2 

Gn=— f,  (3.2) 


where 


G  =  -jtfV2  M  +  j  y/Jiu  V  D  +  K, 


(3.3) 
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and 


M= 


rrih/m  x, 


xn 


K  = 


(Uh/uJa)2  0 
0  rl 


'a  1  rw  _  r (h(uh/Ua)  0  ] 

*a\'  1  0 

•“-RV-te}- 


(3.4) 


3.2.3  Static  Aeroelastic  Equation 

For  the  present  investigation,  since  nonsymmetric  airfoil  sections  and/or  nonzero  angles-of- 
attack  are  being  considered,  one  must  also  include  the  effects  of  changes  in  the  computed 
mean  flow  (e.g.  zeroth  harmonic)  for  different  LCO  amplitudes,  which  in  turn  effect  the 
mean  angle-of-attack  a0  of  the  LCO.  This  requires  that  one  must  also  satisfy  the  static  (e.g. 
zeroth  harmonic)  aeroelastic  equation  for  the  pitch  degree-of-freedom 

Ka(oi0  -  aeo)  =  qoc^scmoiu,  a0,  «i,  h/b),  (3.5) 


where  aeo  is  the  angle-of-attack  for  zero  spring  stiffness  in  the  pitch  coordinate.  After  some 
rearrangement,  this  equation  can  be  written  as 

(3.6) 

The  static  aeroelastic  equation  for  the  plunge  coordinate  ho  establishes  the  vertical  position 
of  the  airfoil.  This  equation  is  decoupled  from  both  the  static  pitch  equation  and  the  first 
harmonic  unsteady  aeroelastic  equations  (Eq.  3.2)  and  thus  is  neglected  in  the  following 
theoretical  development. 


2V2_ 

^eo  o  ^hno  • 

tttI 


3.2.4  LCO  Solution  Procedure 

In  the  following,  we  describe  the  technique  in  which  the  harmonic  balance  method  can  be 
used  for  modeling  limit  cycle  oscillations.  The  methodology  was  initially  developed  and  pre¬ 
sented  in  Thomas  et.  al.  [45].  The  HB/LCO  solution  procedure  starts  by  rewriting  Eq.  3.2  as 

V2 

Gv  =  — —  f,  (3.7) 

7TG!i  ' 

where 
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In  this  form,  Eq.  3.2  has  been  divided  through  by  the  first  harmonic  unsteady  pitch  am¬ 
plitude  ai.  This  enables  one  to  consider  c*i  as  the  independent  variable  in  the  HB/LCO 
solution  process. 

Defining  R(L)  as  the  vector  operator  representing  residual  of  Eq.  3.6  together  with  the 
residual  of  the  real  and  imaginary  parts  of  Eq.  3.7,  one  may  write  the  governing  steady  and 
unsteady  aeroelastic  equations  in  vector  form  as 


R(L)  = 


2V2  _ 

bo  Oie o  „  cmo 

7i TS 
V2  ~ 

Gv  —  —z~f 

7tOii 


>-o. 


(3.9) 


where  G  is  the  4  x  4  matrix 


G  = 


Grte  — Gjm 
Gim  Gr6 


(3.10) 


v  and  f  are 


v  = 


Re(hi/ctib) 
1 

Im(/ii/o:i6) 
0 


f=( 


f  — Re(ci1)  'j 
2Re(cmi) 

[_  2  Im^j )  J 


(3.11) 


and  L  is  the  vector  of  unknown  LCO  variables 


L  =  < 


Oio 

V 

U) 


Re(hi/aib) 

lm(hi/aib)  J 


(3.12) 


We  have  found  that  a  Newton-Raphson  nonlinear  root  finding  technique  is  an  efficient  and 
stable  method  for  quickly  solving  for  a  root  L  of  Eq.  3.9.  That  is,  for  a  specified  unsteady 
pitch  amplitude  <5i,  one  can  implement  an  iterative  process  whereby  the  (n  +  1)^  update 
to  the  LCO  solution  L  is  given  by 


Ln+i  =  Ln 


'<9R(Ln)' 

dL 


R(Ln), 


where 


(3.13) 
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(3.14) 


'dR(L)' 

dL 


9R  &R  9R  8R  <9R 

da0  dV  dui  sRe(/ii/aifc)  8lm(/ii/ai&) 


We  have  also  found  that  one  can  get  very  good  convergence  by  simply  approximating  the 
column  vectors  of  <9R(Ln)  / dL  using  forward  differencing.  That  is, 


dR(L”)  _  R(Ln,  o%  +  e)  -  R(Ln,  6%) 
da0  ~  e 

dR(Ln)  _  R(Ln,  Vn  +  e)-  R(Ln,  Vn) 

■ -  r^Z>  - 


(3.15) 

(3.16) 


etc.  for  a  small  e. 


For  each  step  of  the  LCO  solution  procedure,  the  harmonic  balance  flow  solver  is  imple¬ 
mented  using  the  current  LCO  frequency  u>  and  structural  mode  shape  ( h\/a\b ),  for  the 
prescribed  LCO  pitch  amplitude  di,  to  provide  an  update  for  the  LCO  aerodynamic  load¬ 
ing  (i.e.  q1(  c^,  and  CmJ.  The  technique  is  marched  until  a  suitable  level  of  convergence 
is  achieved,  and  the  linear  flutter  solution  has  been  found  to  provide  an  excellent  starting 
solution  for  the  iterative  process.  Typically,  only  a  few  iterations  are  required  to  achieve 
convergence. 


3.3  Steady  and  Unsteady  Aerodynamic  Modeling  of 
the  NLR  7301  Airfoil 

The  model  configuration  under  consideration  is  the  NLR  7301  constant  airfoil  section  wing 
tested  extensively  by  Schewe  et.  al.  [46,  47]  and  Knipfer  and  Schewe  [48].  Transonic  two- 
degree  of  freedom  aeroelastic  experimental  studies  were  conducted  for  various  Mach  numbers 
and  angles-of-attack,  and  in  some  instances,  LCO  was  observed. 

The  experimental  condition  we  model  using  the  HB/LCO  solution  methodology  is  referred  to 
by  Schewe  et.  al.  [47]  as  ’’measured  point”  condition  MP77.  This  is  an  LCO  condition  where 
the  wind-tunnel  test  section  Mach  number,  mean  angle-of-attack,  and  Reynolds  number  are 
reported  to  be  [48,  47]  —  0.768,  do  =  1-28  degree,  and  Re^  =  1.727  x  106,  respectively. 

The  experimental  model  [46,  48,  47]  consists  of  a  one  meter  span  (s  =  1.0m)  by  0.3  meter 
chord  (c  =  0.3m)  NLR  7301  constant  airfoil  section  wing  placed  in  a  one  meter  by  one  meter 
cross  section  wind  tunnel  test  section.  The  elastic  axis  of  the  model  is  positioned  at  the  wing 
quarter  chord.  Due  to  the  relatively  large  size  of  the  model  in  relation  to  the  wind  tunnel 
cross  section,  wind  tunnel  interference  effects  are  known  to  be  significant. 
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The  harmonic  balance  based  CFD  flow  solver  used  in  the  present  investigation  is  capa¬ 
ble  of  only  modeling  isolated  airfoil  sections  in  an  unbounded  freestream.  As  such,  a  method 
of  accounting  for  wind-tunnel  wall  interference  effects  is  necessary.  Unfortunately,  simple 
classical  analytical  wind-tunnel  wall  interference  correction  methods  have  the  tendency  to 
breakdown  for  transonic  Mach  numbers.  As  such,  one  method  to  account  for  wind-tunnel 
wall  interference  effects  is  to  seek  out  a  combination  of  Mach  number  and  angle-of-attack 
for  the  CFD  method  such  that  the  computed  surface  pressure  distribution  matches  as  best 
as  possible  the  experimental  distribution.  Investigators  such  as  Weber  et.  al.  [54]  and  Tang 
et.  al.  [55],  who  are  also  studying  computational  techniques  for  modeling  test  condition 
MP77,  have  reported  recently  however  that  such  a  combination  of  corrected  Mach  number 
and  angle-of-attack  has  been  quite  elusive  to  determine,  and  appears  perhaps  not  to  exist.  In 
fact,  Castro  et.  al.  [?,  57]  have  recently  taken  the  approach  of  adding  appropriate  boundary 
conditions  to  their  CFD  flow  solver  to  account  for  the  wind-tunnel  walls.  However,  Castro 
et.  al.  have  noted  that  there  is  still  much  difficulty  in  correctly  accounting  for  wind-tunnel 
wall  porosity  effects. 

Rather  than  go  through  a  similar  and  time  consuming  exercise  of  trying  to  determine  a 
combination  of  corrected  Mach  number  and  angle-of-attack  in  an  effort  to  precisely  match 
the  experimental  surface  pressure  distribution,  we  have  decided  as  a  first  approximation  to 
proceed  by  fixing  the  Mach  number  to  =  0.75,  and  to  then  simply  adjust  the  angle-of- 
attack  until  the  computed  CFD  solution  steady  flow  lift  matches  the  experimental  lift. 

3.3.1  Computational  Mesh 

Figure  3.3  shows  sample  inviscid  and  viscous  computational  grids  used  by  the  HB  CFD 
solver  for  the  NLR  7301  airfoil  configuration.  A  ”  C-grid”  structured  mesh  topology  is  used. 
The  ’’fine”  mesh  viscous  grid  shown  in  Fig.  3.3c  and  Fig.  3.3d  consists  of  65  mesh  points 
radially  and  257  mesh  points  circumferentially,  with  a  total  of  193  mesh  points  surrounding 
the  airfoil  surface,  and  the  remainder  of  the  circumferential  mesh  points  being  distributed 
in  the  wake.  For  all  grids,  the  outer  boundary  extends  to  a  distance  of  10  chord  lengths 
from  the  center  of  the  airfoil.  Similar  129  x  33  and  193  x  49  meshes  have  also  been  created. 
Figure  3.3a  and  Fig.  3.3b  show  an  inviscid  193  x  49  mesh. 

For  viscous  simulations,  the  grid  spacing  has  been  placed  close  enough  to  the  airfoil  sur¬ 
face  such  that  the  maximum  computed  value  for  the  nondimensional  distance  to  the  wall 
spacing  parameter  y+  is  everywhere  less  than  2.0  for  the  three  mesh  resolutions  considered. 

3.3.2  Steady  Flow  Simulations 
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LCO  Amplitude  Reduced  Velocity,  V 


Freestream  Mach  Number,  Mto 


(a)  Mach  Number  Flutter  Speed  Dip  Trend 
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(b)  LCO  Behavior  Trends  Near  Flutter  Onset  Condition 

Figure  3.1:  Example  LCO  Behavior  Trends. 
42 


Figure  3.2:  Geometry  for  ”  Typical”  (Pitch/Plunge)  Two-Degree-of- Freedom  Airfoil  Section 
Aeroelastic  Model. 
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(a)  Inviscid  Mesh  Close-Up 


(b)  Inviscid  Mesh  Overall 


(c)  Viscous  Mesh  Close-Up 


Figure  3.3:  Computational  Grids  Used  for  the  NLR  7301  Airfoil  Configuration. 
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Figure  3.4:  Computed  Steady  Flow  Surface  Pressure  Distributions:  NLR  7301  Airfoil  Sec 
tion,  Mqo  =  0.75,  Re oo  =  1.727  X  106. 


Number  of  Harmonics 

Mesh  Size 

1 

2 

3 

129  x  33 
193  x  49 
257  x  65 

(-0.18,-0.40) 

(-0.11,-0.54) 

(-0.08,-0.60) 

(-0.31,-0.29) 

(-0.34,-0.44) 

(-0.32,-0.52) 

(-0.27,-0.26) 

(-0.28,-0.39) 

(-0.28,-0.47) 

Table  3.1:  Mesh  Resolution  and  Number  of  Harmonic  Used  in  Harmonic  Balance  Solver 
Trends  for  Computed  First  Harmonic  Unsteady  Moment  Coefficient  cmi/6i\:  =  0.75, 

Re oo  =  1.727  x  106,  do  =  0.2  (deg),  di  =  5.  (deg),  u  =  0.3,  and  a  =  —0.5. 


With  the  Mach  number  set  to  =  0.75,  steady  flows  are  computed  by  simply  running 
the  HB  flow  solver  with  zero  harmonics.  In  order  to  match  the  experimentally  observed  lift 
coefficient  of  c/0  =  0.27,  an  angle-of-attack  of  do  =  0.2  degree  is  needed  for  the  viscous  flow 
model,  while  an  angle-of-attack  of  do  =  —1.55  degree  is  required  for  the  inviscid  model. 

For  both  the  inviscid  and  viscous  flow  models,  Fig.  3.4  shows  the  computed  and  experi¬ 
mental  steady  surface  pressure  distributions,  while  Fig.  3.5  shows  computed  steady  Mach 
contours.  In  this  instance,  a  257  x  65  grid  has  been  used  for  the  viscous  model,  and  a  193  x  49 
grid  has  be  used  for  the  inviscid  model.  A  strong  shock  is  readily  apparent  for  the  inviscid 
model. 


3.3.3  Unsteady  Flow  Simulations 

Next,  Fig.  3.6  shows  computed  Mach  contours  during  one  cycle  of  motion  for  the  fine  mesh 
(257  x  65)  viscous  flow  CFD  model  of  the  NLR  7301  airfoil  when  oscillating  about  the 
quarter-chord  at  a  pitch  amplitude  of  <5i  =  5.  degrees  and  a  reduced  frequency  u>  =  0.3. 
Clearly  evident  is  shock  induced  boundary  layer  separation  over  both  the  upper  and  lower 
surfaces  during  a  cycle  of  motion. 

3.3.4  Mesh  and  Harmonic  Convergence  Issues 

To  have  confidence  in  the  LCO  solutions,  one  needs  to  determine  whether  or  not  a  sufficient 
level  of  grid  resolution  is  being  used  to  model  the  correct  flow  physics  accurately.  For  the 
harmonic  balance  solver,  the  issue  of  solution  resolution  is  further  complicated  by  the  fact 
that  one  must  also  consider  harmonic  convergence.  That  is,  whether  or  not  a  sufficient  num¬ 
ber  of  harmonic  expansion  terms  iV#  are  being  used  in  the  HB  flow  solver. 

Thus,  much  like  when  working  with  time  domain  methods  where  one  must  also  be  con¬ 
cerned  with  temporal  accuracy  in  addition  to  spatial  accuracy,  for  the  HB  method,  both 
mesh  and  harmonic  resolution  must  be  considered.  To  help  give  a  sense  if  sufficient  mesh 
resolution  and  harmonics  are  being  utilized,  one  can  define  a  ’’most  difficult”  test  case.  For 
example,  consider  a  configuration  that  is  moving  at  the  largest  amplitude  that  one  wishes 
to  model,  and  also  near  the  reduced  frequency  that  is  of  interest,  say  the  flutter  reduced 
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Dimensional  Quantities 

m 

=  4.65  kg 

mSD 

=  8.24  kg 

mSp 

=  13.75  kg 

la 

=  0.086  kg  m? 

Kh 

=  1.21  x  106  N/m 

Ka 

—  6.68  x  103  N  m/rad 

Sa 

=  0.387  kg  m 

Dh 

=  82.9  kg /sec 

Da 

=  0.197  kg  m2/ (sec  rad ) 

Poo 

=  0.383  kg/m 3 

Nondimensional  Quantities 

V 

=  172 

•E a. 

=  0.555 

r 2 
'  a 

=  0.822 

Uh/Ua 

=  1.83 

Ch 

=  0.0175 

(a 

=  0.00411 

Table  3.2:  Structural  Parameter  Values  for  the  NLR  7301  Aeroelastic  Configuration. 

frequency.  One  can  then  create  a  table  of  computed  unsteady  aerodynamic  lift  or  moment 
values  for  different  mesh  resolutions  and  numbers  of  harmonics  (NH)  used  in  the  HB  flow 
solver. 

Table  3.1  in  fact  shows  a  total  of  nine  individually  computed  values  for  unsteady  moment 
for  the  three  different  mesh  resolutions  (129  x  33,  193  x  49,  257  x  65)  and  either  one,  two,  or 
three  harmonics  (iV#)  used  in  the  HB  algorithm.  An  unsteady  pitch  amplitude  of  di  =  5. 
degrees  is  chosen,  which  is  the  maximum  amplitude  we  consider  in  the  LCO  analysis  to  fol¬ 
low,  along  with  a  reduced  frequency  of  u>=0.3.  As  can  be  seen,  the  unsteady  solution  for  this 
’’most  difficult”  test  case  appears  to  be  rapidly  converging.  Therefore,  for  the  LCO  studies 
presented  in  the  following,  results  are  based  on  the  193  x  49  mesh  using  two  harmonics 
(Nh  =  2)  in  the  HB  flow  solver. 

3.4  Limit  Cycle  Oscillation  Modeling  for  the  NLR  7301 
Configuration 

3.4.1  Structural  Parameters 

The  structural  parameters  used  for  the  LCO  analysis  are  those  presented  in  Knipfer  and 
Schewe  [48]  and  Schewe  et.  al.  [47].  Table  3.2  summarizes  values  for  dimensional  and 
nondimensional  quantities. 
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3.4.2  Flutter  Point  Prediction 


Computing  the  LCO  behavior  via  the  HB/LCO  solution  methodology  proceeds  by  first  deter¬ 
mining  the  flutter  onset  condition.  As  noted  before,  experimental  measured  point  condition 
MP77  corresponds  to  an  LCO  condition.  However,  the  observed  LCO  amplitude  is  very 
small.  Approximately  |dx|  =  0.18  degree  (Fig.  15  of  Knipfer  and  Schewe  [48]).  As  will  be 
shown  subsequently,  based  on  the  predicted  LCO  behavior  using  the  HB/LCO  method,  for 
such  a  small  amplitude,  the  flow  conditions  at  measured  point  MP77  are  likely  to  be  very 
near  the  flutter  onset  condition.  For  this  analysis,  we  are  considering  this  to  be  the  case.  As 
such,  we  are  in  effect  assuming  that  the  lift  coefficient  at  the  flutter  onset  condition  is  also 
ci0f  =  0.27. 

Rather  than  use  a  separate  linearized  unsteady  aerodynamic  solver  to  determine  the  flutter 
onset  condition,  one  can  instead  simply  run  the  HB  CFD  flow  solver  with  a  single  harmonic 
(Nh  =1)  and  consider  a  very  small  amplitude  motion.  As  the  amplitude  of  motion  goes  to 
zero,  the  HB  solver  yields  exactly  the  same  answer  as  a  linearized  unsteady  solver.  This  is 
another  useful  feature  of  the  HB  solver. 

Predicting  the  linear  flutter  point  condition  thus  consists  of  first  tabulating  linear  unsteady 
aerodynamic  loading  data  for  plunge  and  pitch  motions,  respectively,  for  discrete  reduced 
frequencies  cD  over  a  range  of  reduced  frequencies  where  flutter  is  presumed  to  occur.  In  this 
instance,  we  have  chosen  small  unsteady  motion  amplitudes  of  h\/b  =  0.001  and  ot\  =  0.001 
degree,  respectively  for  the  plunge  and  pitch  degrees-of-freedom,  in  order  to  simulate  dynam¬ 
ically  linear  unsteady  aerodynamics.  Normalizing  the  calculated  unsteady  lift  and  moment 
coefficients  by  the  respective  amplitudes  of  either  the  plunge  or  the  pitching  motions,  one 
then  can  obtain  tabulated  data  for  the  unsteady  aerodynamic  transfer  functions  Q1Ri/6(u>), 
q1_i  (Q),  cmi-i/t  (u>),  and  cmi_i  (d>)  at  the  preselected  reduced  frequencies  u).  The  2x2  linear 
unsteady  aerodynamic  transfer  function  matrix  F(u>)  is 

F  =  [  C‘iSl(“)  1 

|_  ^1/^/6  M  CnUaiM]  ’ 

and  for  dynamically  linear  unsteady  aerodynamics, 

f  =  Fu.  (3.18) 

With  the  linear  unsteady  aerodynamic  transfer  functions  tabulated  for  a  range  of  real  re¬ 
duced  frequencies,  one  can  then  use  a  ”  V-g”  or  similar  type  of  method  to  solve  for  the  flutter 
onset  conditions.  For  example,  by  rewriting  Eq.  3.2  as 

(3.19) 

one  can  search  the  plane  of  values  for  reduced  velocity  and  reduced  frequency,  (V,Q),  for  the 


f  G  -  u  =  H(V,  d>)u  =  0, 


(3.17) 
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Variable 

Inviscid 

Viscous 

Experiment 

50  (deg) 

-1.55 

0.20 

Qo 

0.267 

0.272 

0.27 

C77l0 

-0.140 

-0.0646 

-0.082 

V 

0.371 

0.385 

0.488 

U) 

0.312 

0.301 

0.242 

Oil  (deg) 

0.0001 

0.0001 

0.18 

{hx/otib) 

(1.15,0.141) 

(1.12, .0.172) 

(1.66,0.) 

Table  3.3:  Computed  Flutter  Point  Conditions  and  NLR  7301  Aeroelastic  Configuration 
Experimental  Test  Case  MP77  LCO  Conditions. 


condition(s)  where  the  magnitude  of  the  determinant  of  the  2x2  matrix  H  is  a  minimum. 
Using  this  technique  provides  a  quick  and  easy  way  of  determining  a  good  approximation 
of  the  flutter  onset  condition(s).  This  approximate  flutter  solution  then  provides  a  good 
starting  solution  for  the  HB/LCO  method,  which  when  then  run  for  a  very  small  amplitude 
motion,  can  be  used  to  seek  out  the  precise  flutter  onset  condition.  This  is  what  has  been 
done  for  the  current  flutter  analysis  of  the  NLR  7301  configuration,  and  the  results  for  the 
precise  flutter  onset  conditions  for  the  viscous  and  inviscid  flow  models  are  shown  in  Ta¬ 
ble  3.3,  which  also  presents  data  for  the  LCO  condition  of  experimental  test  point  MP77. 

As  can  be  seen  from  Table  3.3,  the  computed  flutter  onset  conditions  are  somewhat  different 
than  those  observed  for  the  MP77  test  case.  The  numerical  model  predicts  a  lower  flutter 
reduced  velocity,  and  a  higher  flutter  reduced  frequency,  respectively  for  both  the  inviscid 
and  viscous  flow  models.  The  MP77  experimental  result  for  the  structural  eigenvector  is 
also  somewhat  more  plunge  dominated  than  the  computed  values.  Figure  15  of  Knipfer 
and  Schewe  [48]  was  used  to  determine  the  value  of  (hi/atib)  for  the  MP77  experimental 
condition.  Note,  Knipfer  and  Schewe  [48]  define  the  plunge  coordinate  h  positive  upward, 
opposite  to  our  definition.  We  have  corrected  for  this  sign  difference  for  the  experimental 
value  listed  in  Table  3.3. 

We  are  currently  trying  to  ascertain  what  possible  causes  might  explain  the  differences 
between  the  experimental  and  computed  conditions.  Perhaps  the  steady  flow  CFD  solu¬ 
tion  based  on  adjusting  the  angle-of-attack  to  match  the  experimental  lift  coefficient  does 
not  provide  a  sufficient  approximation  of  the  true  wind-tunnel  aerodynamic  environment,  or 
perhaps  the  MP77  test  case  corresponds  to  a  flow  condition  more  removed  from  the  flutter 
onset  condition. 

We  wish  to  also  point  out  that  we  define  mass  ratio  //  based  on  the  wing  model  mass 
m,  whereas  Schewe  et.  al.  [46,  47]  and  Knipfer  and  Schewe  [48]  define  the  mass  ratio  based 
on  the  total  plunge  coordinate  mass  m^.  This  in  turn  means  that  the  reduced  velocity  is 
also  defined  somewhat  differently.  One  may  relate  the  two  definitions  for  reduced  velocity 
via  ^Jirrih/m).  That  is,  our  definition  for  reduced  velocity  is  the  product  of  the  Schewe  et. 
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al.  [46,  47]  and  Knipfer  and  Schewe  [48]  definition  multiplied  by  y^m^/ra). 

3.4.3  Zero  Spring  Stiffness  Angle-of- Attack 

Next,  the  angle-of-attack  corresponding  to  zero  torsional  spring  stiffness,  aeo,  is  computed 
for  both  the  viscous  and  inviscid  models.  Once  the  flutter  point  is  established,  this  value 
is  easily  determined  using  Eq.  3.6.  For  the  inviscid  flow  model,  the  value  is  aeo  =  —0.692 
(deg),  and  for  the  viscous  flow  model,  the  value  is  aeo  =  0.645  (deg). 

3.4.4  Computed  LCO  Behavior 

With  the  flutter  onset  condition  established,  along  with  a  computed  value  for  the  angle-of- 
attack  for  zero  torsional  spring  stiffness  aeo,  one  can  proceed  with  the  LCO  analysis  for  finite 
amplitude  motions.  Figure  3.7  shows  computed  LCO  solution  behavior  trends  for  both  the 
viscous  and  inviscid  models  of  the  NLR  7301  configuration.  Also  shown  in  Figure  3.7  are  the 
experimental  LCO  conditions  as  reported  by  Schewe  et.  al.  [46,  47]  along  with  computed 
LCO  conditions  as  determined  by  the  investigations  of  Weber  et.  al.  [54]  and  Tang  et. 
al.  [55]. 

Clearly  from  the  LCO  amplitude  versus  reduced  velocity  curves  shown  in  Fig.  3.7a,  modeling 
viscous  flow  effects  leads  to  a  nonlinear  soft  or  benign  LCO  behavior  trend  as  discussed  in  the 
introduction,  whereas  the  inviscid  model  exhibits  a  nearly  linear  aeroelastic  LCO  behavior 
trend.  The  flutter  onset  condition  is  also  indicated  in  the  figures  and  is  less  sensitive  to 
viscous  effects.  Fig.  3.7b  and  Fig.  3.7c  show,  respectively,  LCO  behavior  trends  for  the 
frequency  ratio  and  reduced  frequency.  The  reduced  frequency  can  be  seen  to  change  by 
a  significant  amount  in  the  viscous  case  for  large  amplitudes.  Finally,  Fig.  3.7d  shows  the 
change  in  the  LCO  mean  angle-of-attack  do  with  respect  to  reduced  velocity.  The  LCO 
mean  angle-of-attack  is  observed  to  decrease  by  approximately  0.5  degree  from  the  flutter 
onset  condition  for  the  range  of  LCO  amplitudes  considered  in  the  case  of  the  viscous  model. 
Since  the  aerodynamics  are  known  to  be  extremely  sensitive  to  small  changes  in  the  transonic 
region,  it  is  likely  to  be  important  to  model  this  effect. 
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(b)  to  +  AT 


(c)  to  +  2  AT 


(d)  t0  +  SAT 


(e)  to  +  4AT 


(g)  to  +  6A T 


Figure  3.6:  Computed  HB  Method  Unsteady  Mach  Number  Contours:  NLR  7301  Airfoil 
Section,  =  0.75,  Re oo  =  1.727  X  106,  6lq  =  0.2  (deg),  «i  =  5.  (deg),  a>  =  0.3, 
a  =  -0.5,  Nh  =  3  (AT  =  T/(2NH  +  1  f%=  T/7). 
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Figure  3.7:  LCO  Behavior  Trends  for  the  NLR  7301  Configuration 
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3.5  Conclusions 


An  efficient  procedure  for  computing  the  LCO  behavior  of  aeroelastic  airfoil  configurations 
in  viscous  transonic  flows  is  presented.  Viscous  effects  are  demonstrated  to  be  an  important 
factor  in  determining  the  behavior  of  LCO  response  with  respect  to  reduced  velocity.  Vis¬ 
cous  effects  for  the  NLR  7301  aeroelastic  configuration  under  investigation  lead  to  a  much 
more  gradual  rate  of  increase  in  LCO  amplitude  with  respect  to  reduced  velocity  when  com¬ 
pared  to  an  inviscid  model.  The  rapid  increase  in  the  computed  LCO  response  beyond  the 
flutter  point  emphasizes  the  importance  of  making  careful  and  comprehensive  experimental 
measurements  of  LCO  response  over  a  range  of  reduced  velocities  or  other  parameters  as  the 
flutter  boundary  is  exceeded. 


Nomenclature 


a 

6,  C 


Cfo  )  Cjtjq 


ch  >  Cjtii 


Dh,Da 


e 

h:  a 


ho ,  olq 


h-i.a-i 


Ia 

3 

Kh,Ka 

Moo 

m 


=  nondimensional  location  of  airfoil 
elastic  axis,  a  =  e/6  (see  Fig.  3.2) 

=  semi-chord  and  chord,  respectively 
=  coefficient  of  lift  and  moment  about 
elastic  axis,  respectively 
=  zeroth  harmonic  or  mean  coefficient 
of  lift  and  moment  about 
elastic  axis,  respectively 
=  first  harmonic  unsteady  coefficient 
of  lift  and  moment  about 
elastic  axis,  respectively 
=  airfoil  plunge  damping  and  torsional 
damping,  respectively 
=  location  of  airfoil  elastic  axis,  measured 
positive  aft  of  airfoil  midchord 
=  airfoil  plunge  and  pitch  coordinate 
degree  of  freedom,  respectively 
=  zeroth  harmonic  or  mean  airfoil  plunge 
and  pitch  amplitude,  respectively 
=  first  harmonic  unsteady  airfoil  plunge 
and  pitch  amplitude,  respectively 
=  second  moment  of  inertia  of  airfoil 
about  elastic  axis 

=  V'-T 

=  airfoil  plunge  stiffness  and  torsional 
stiffness  about  elastic  axis,  respectively 
=  freestream  Mach  number 
=  wind-tunnel  experimental  wing  mass 
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m-SD  =  wind-tunnel  shift  device  mass 
msp  =  wind-tunnel  suspension  device  mass 
nth  —  wind-tunnel  plunge  coordinate  total 
mass,  mh  =  m  +  msD  +  rnSP 
p  =  mass  ratio,  p  =  m/'KpO0b2s 

v  =  kinematic  viscosity 

Nh  =  number  of  harmonics  used  in  harmonic 
balance  CFD  flow  solver  method 
Re oo  =  freestream  Reynolds  number 

ra  =  radius  of  gyration  of  airfoil  about 
elastic  axis,  r2  =  Ia/mb2 

Pi  Poo  —  local  and  freestream  density,  respectively 
Sa  =  first  moment  of  inertia  of  airfoil  about 
elastic  axis 

s  =  wind-tunnel  experimental  model  span 
T  =  period,  time  for  one  cycle  of  oscillation 
tw  =  airfoil  surface  tangential  shear  stress 
T  =  period,  time  for  one  cycle  of  oscillation 
Uoo  =  freestream  velocity 

V  =  reduced  velocity,  V  =  U^/y/p, uab 
u,u)  =  frequency  and  reduced  frequency  based 
on  airfoil  chord  Q  =  o;c/C700,  respectively 
=  plunge  coordinate  natural  frequency 
based  on  wing  mass,  uh  =  y/Kh/m 
ua  =  pitch  coordinate  natural  frequency, 
wa  =  \/Ka/Ia 

xa  =  airfoil  static  unbalance,  xa  =  Sa/mb 
y  =  dimensional  distance  normal  to  the 

airfoil  surface 

y+  =  nondimensional  distance  normal  to  the 
airfoil  surface,  y+  =  yy/xjp/u 
Ch  =  plunge  coordinate  damping  coefficient, 

Ch  =  Dh/2muh 

Ca  =  pitch  coordinate  damping  coefficient, 

Co: 


Subscripts/Superscript 

/  =  flutter  point  condition 

Re,  Im  =  real  and  imaginary  part,  respectively 
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Chapter  4 


ZONA  2-DOF  Airfoil  LCO: 
CFL3D.AE 


Summary 

•  Several  efforts  were  directed  toward  the  validation  of  the  measured  LCO/flutter  data 
of  DLR  and  the  verification  of  the  results  due  to  Duke  University  using  the  frequency- 
domain  POD/ROM  eigenmode  approach.  In  Phase  II,  the  work  in  this  chapter  using 
a  time-marching  CFD  method  in  two  dimensions  for  the  simulation  of  the  transonic 
LCO/flutter  of  a  supercritical  airfoil  (NLR  7301)  are  presented  in  three  parts 

1.  Fundamentals  and  preliminary  formulations  and  result  comparison  with  DLR’s 
measured  data  which  is  presented  in  this  chapter. 

2.  Further  refinement  of  (1)  in  various  details,  which  is  published  in  the  Journal  of 
Fluid  and  Structures  (17  (2003)29-41 )■ 

3.  Extension  of  parts  (1)  and  (2)  in  the  investigation  on  the  impact  of  initial  condi¬ 
tion  on  the  transonic  LCO  to  the  same  NLR  airfoil,  which  is  presented  in  Appendix 
A. 

The  work  in  this  part  is  necessary,  as  such  impact  on  LCO  by  initial  conditions  could 
not  be  so  investigated  by  the  frequency- domain  approach  of  Duke  using  harmonic  bal¬ 
ance  (HB)  technique.  ( This  work  is  presented  in  Symposium  Transsonicum  IV  held  in 
Goettingen,  Germany,  Sept.  2002  and  is  published  as  Symposium  Transsonicum  IV  by 
Kluwer  Academic  Publishers  Nov.  2003  (ISBN1-4020-1608-5)) 

•  CFD-based  aeroelastic  computation  is  performed  to  investigate  the  effect  of  nonlin¬ 
ear  aerodynamics  on  transonic  LCO  (Limit  Cycle  Oscillations)  characteristics  of  a 
two-dimensional  supercritical  wing  with  the  NLR  7301section.  It  is  found  that  LCO 
amplitude  given  by  Navier-Stokes  computation  is  less  than  1/3  of  the  one  given  by  Eu¬ 
ler  computation,  much  closer  to  the  experimental  data.  But  Navier-Stokes  computation 
of  LCO  is  turbulence  model  dependent.  Baldwin-Lomax  model  with  Degani-Schiff  mod¬ 
ification  seems  to  produce  spurious  vorticity  for  moving  grids.  Contrary  to  common 
sense,  a  larger  amplitude  initial  perturbation  produces  an  smaller  amplitude  LCO  for 
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Figure  4.1:  Transonic  LCO  Solutions  and  DLR  Test  Data  for  the  NLR7301  Supercritical 
Airfoil  at  M=0.75. 

the  present  case.  The  mean  positions  may  change  with  the  amplitude  of  the  initial  per¬ 
turbation.  Also  addressed  in  the  paper  are  certain  aspects  of  multiblock  MPI  (Message 
Passing  Interface)  parallel  computation  related  to  the  present  problem.  Taken  from 
the  abstract  of  the  paper  entitled  ’’Numerical  Investigation  of  Transonic  Limit  Cycle 
Oscillations  of  a  Two  Dimensional  Supercritical  Wing.”  (above  part  (2)) 

•  A  CFD  time-marching  method,  CFL3D  v6,  is  employed  to  investigate  the  influences 
of  viscosity,  initial  airfoil  position,  and  initial  perturbation  on  the  amplitudes  of  the 
transonic  limit  cycle  oscillations  (LCO)  of  a  supercritical  airfoil.  As  expected,  stronger 
aerodynamic  nonlinearity  leads  to  smaller  LCO  amplitudes,  even  a  damped  solution 
while  weaker  aerodynamic  nonlinearity  incurs  larger  LCO  amplitudes,  even  a  divergent 
solution.  Taken  from  the  abstract  of  the  paper  entitled,  ’’Nonlinear  Aerodynamic  Effects 
on  Transonic  LCO  Amplitude  of  a  Supercritical  Airfoil.”  (above  part  (3)) 

•  The  computed  results  due  to  different  methods  on  the  transonic  LCO  solution  for 
NLR7301  supercritical  airfoil  at  M  =  0.75  can  be  summarized  by  the  following  fig¬ 
ure  prepared  by  Duke  University.  It  is  seen  that  the  DLR  measured  LCO  amplitude 
(Schewe)  is  only  a  fraction  of  the  computed  results  including  that  of  ZONA  (Tang  et 
al),  Naval  Postgraduate  School  (Weber  et  al)  and  Duke  University  (Hall  and  Dowell, 
Hollow  Symbols). 


4.1  Introduction 

LCO  has  been  a  persistent  problem  on  several  current  fighter  aircraft  and  is  generally  en¬ 
countered  with  external  store  configurations.  Denegri  (2000)  provided  a  detailed  description 
of  the  aircraft/store  LCO  phenomenon.  Norton  (1990)  gave  an  excellent  overview  of  LCO 
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of  fighter  aircraft  carrying  external  stores  and  its  sensitivity  to  store  carriage  configuration 
and  mass  properties. 

LCO  can  be  characterized  as  sustained  periodic  oscillations  which  neither  increase  nor  de¬ 
crease  in  amplitude  over  time  for  a  given  flight  condition.  Using  an  s-domain  unsteady  aero¬ 
dynamic  model  of  the  aircraft  and  stores,  Chen  et  al.  (1998)  have  shown  that  wing/store 
LCO  can  be  a  post-flutter  phenomenon  whenever  the  flutter  mode  contains  low  unstable 
damping.  This  type  of  flutter  mode  is  called  a  ’’hump  mode”.  Since  the  aircraft  structure 
usually  contains  structural  nonlinearity  such  as  friction  damping,  this  amplitude-dependent 
friction  damping  can  suppress  the  growth  of  amplitude,  thus  resulting  in  a  steady  state  os¬ 
cillation.  This  is  known  as  the  nonlinear  structural  damping  (NSD)  model  of  the  wing/store 
LCO.  Although  not  thoroughly  proven  through  tests  or  numerical  simulations,  results  of  the 
NSD  show  excellent  correlation  with  flight  test  LCO  data  of  F-16  throughout  subsonic  and 
transonic  Mach  numbers.  On  the  other  hand,  other  researchers,  notably  Cunningham  and 
Meijer  (1995),  believe  that  the  wing/store  LCO  is  due  largely  to  the  transonic  shock  oscil¬ 
lation  and  shock  induced  flow  separation,  called  Transonic  Shock/Separation  (TSS)  model. 
Edwards  has  suggested  the  TSS  model  and  viscous  effects  are  two  major  factors  that  cause 
transonic  LCO  for  wings.  He  also  has  studied  the  shock  buffet  phenomenon  in  addition  to 
transonic  LCO  [65].  It  should  be  noted  that,  however,  there  is  no  conflict  in  the  NSD  model 
and  the  TSS  model  in  that  both  physical  effects  may  contribute  to  the  various  instances  of 
LCO. 

Recent  renewed  interest  in  LCO  is  perhaps  motivated  by  the  need  to  further  understand  the 
physics  of  LCO  and  the  current  advent  of  CFD  methodology  in  aeroelasticity.  Among  the 
potential  computational  methods  for  LCO  prediction/investigation,  two  will  be  mentioned 
here:  the  CFL3D  code  (version  6)  [59],  [69,  72],  developed  and  supported  by  NASA/Langley 
and  the  POD/ROM  EigenMode  approach  [67],  originated  by  Dowell  and  Hall  of  Duke  Univer¬ 
sity.  The  former  is  a  conventional  time-domain  CFD  method  whereas  the  latter  a  frequency- 
domain  CFD  method,  using  aerodynamic  eigenmodes. 

The  present  study  uses  the  CFD  time-marching  method,  CFL3D  v6,  to  numerically  inves¬ 
tigate  transonic  LCO  of  a  two-dimensional  supercritical  wing  under  a  plunging/ pitching 
spring-mounting  system  [70,  71,  75].  It  is  reasonable  to  start  by  investigating  a  two- 
dimensional  LCO  case  in  order  to  better  understand  the  physics  of  LCO.  However,  because 
of  the  complexity  of  a  two-dimensional  LCO  experimental  test,  there  are  few  experimental 
data  available  for  comparison.  To  the  best  of  our  knowledge,  the  experimental  work  per¬ 
formed  by  Schewe  et.  al.  [70,  71,  75],  is  perhaps  the  only  two-dimensional  LCO  experimental 
test  available  in  the  literature.  Those  test  data  were  immediately  used  by  Platzer  et.  al. 
to  validate  their  thin-layer  Navier-Stokes  aeroelastic  solver  [60,  77].  While  the  emphasis  of 
[60,  77],  was  on  the  predictive  capability  of  the  thin-layer  Navier-Stokes  aeroelastic  solver, 
our  emphasis  here  is  to  investigate  the  effect  of  aerodynamic  viscosity,  turbulence  modeling 
and  solution  strategies,  such  as  initial  perturbation  size,  on  transonic  LCO  of  the  supercrit¬ 
ical  airfoil.  Also  addressed  in  the  paper  are  some  of  the  issues  related  to  multiblock  MPI 
parallel  aeroelastic  computation  of  this  problem. 
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4.2  Numerical  Methodology 

The  computer  code  used  in  this  study  is  CFL3D  v6,  which  solves  the  three-dimensional 
thin-layer  Reynolds  averaged  Navier-Stokes  equations  with  an  upwind  finite  volume  formu¬ 
lation  [72].  A  two-dimensional  problem  can  be  calculated  by  using  two  identical  grid  planes, 
created  by  duplicating  the  two-dimensional  grid. 

The  code  uses  formally  third-order  upwind-biased  spatial  differencing  for  the  inviscid  terms 
with  flux  limiting  in  the  presence  of  shocks.  Either  flux-difference  or  flux-vector  splitting 
are  available.  The  flux-difference  splitting  method  of  Roe  (1981)  is  employed  in  the  present 
computations  to  obtain  fluxes  at  cell  faces.  Viscous  terms  are  discretized  with  second-order 
central  differencing.  There  are  two  types  of  time  discretization  available  in  the  code.  The 
first-order  backward  time  differencing  is  used  for  steady  calculation  while  the  second-order 
backward  time  differencing  with  t-TS  subiterations  is  used  for  static  and  dynamic  aeroelastic 
calculation.  Furthermore,  grid  sequencing  for  steady  state  and  multigrid  and  local  pseudo¬ 
time  stepping  for  time  marching  solutions  are  employed.  Also  available  in  the  code  are  many 
turbulence  models,  although  here  only  the  Spalart-Allmaras  model  [76],  and  Baldwin-Lomax 
model  [58],  with  the  Degani-Sehiff  modification  have  been  used.  A  detailed  description  of 
the  methodology  of  the  code  can  be  found  in  [72]. 

One  of  the  important  features  of  the  CFL3D  code  is  its  capability  of  solving  multiple  zone 
grids  with  one-to-one  connectivity.  Spatial  accuracy  is  maintained  at  zone  boundaries,  al¬ 
though  subiterative  updating  of  boundary  information  is  required.  Coarse-grained  paral¬ 
lelization  using  the  MPI  protocol  can  be  utilized  in  multiblock  computations  by  solving  one 
or  more  blocks  per  processor.  When  there  are  more  blocks  than  processors,  optimal  perfor¬ 
mance  is  achieved  by  allocating  an  equal  number  of  blocks  to  each  processor.  As  a  result, 
the  time  required  for  a  CFD-based  aeroelastic  computation  can  be  dramatically  reduced. 
In  this  paper,  both  single  and  multiblock  MPI  parallel  aeroelastic  computations  near  the 
onset  of  flutter  LCO  are  compared  with  experiment  and  with  other  computations.  Figure 
4.2  shows  a  C-type  grid  with  273  x  93  mesh  points  around  the  NLR  7301  airfoil  that  has 
been  divided  into  eight  69x47  blocks.  This  and  a  single  block  version  of  this  grid  are  used 
in  the  computations  to  follow. 

In  order  to  accommodate  multiblock  computation,  the  mesh  deformation  scheme  in  [59] 
is  modified.  In  [59],  the  remeshing  scheme  uses  a  modified  spring  analogy  with  solid  body 
translation/rotation  of  the  fluid  mesh  near  solid  surfaces.  Initialization  of  the  grid  defor¬ 
mation  at  each  step  is  performed  using  a  TFI  (Transfinite  Interpolation)  step.  The  mesh 
interior  is  then  smoothed  and  grid  orientation  near  boundaries  is  preserved  using  the  modi¬ 
fied  spring  analogy. 

In  the  present  implementation,  the  subgrid  based  TFI  scheme  of  [68],  has  been  employed  for 
initialization  at  each  time  step.  That  scheme  uses  subgrids  consisting  of  ’’slave  vertices”  to 
move  both  block  boundaries  and  interiors.  In  some  instances,  in  order  to  achieve  an  optimal 
division  of  grid  points,  it  is  necessary  to  place  flow  field  block  boundaries  near  a  moving  solid 
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Figure  4.2:  Multiblock  C-type  grid  around  NLR  7301  airfoil  (eight  69x47  blocks). 


surface.  An  example  of  this  is  shown  in  Figure  4.2.  The  multiblock  boundary  and  interior 
movement  scheme  allows  the  user  to  place  block  boundaries  near  surfaces  as  necessary  for 
optimal  parallelization.  Boundaries  interior  to  the  fluid  domain  near  a  surface  respond  to 
the  local  surface  motion.  As  the  airfoil  moves,  block  boundaries  move  to  maintain  integrity 
of  block  interfaces  and  the  airfoil  surface. 

User  controlled  input  makes  it  possible  to  update  the  mesh  using  this  subgrid/TFI-based 
scheme  alone  or  to  update  with  an  initialization  using  this  scheme  plus  additional  smoothing 
steps.  These  added  smoothing  steps,  the  number  of  which  can  be  defined  by  the  user,  employ 
the  modified  spring  analogy  scheme  [59].  In  the  current  implementation  the  spring  analogy 
scheme  updating  the  mesh  interior  is  now  written  in  delta  formulation  so  that  the  relative 
orientation  of  the  original  grid  is  retained.  The  solid  body  rotation/translation  of  the  fluid 
grid  is  also  now  performed  near  both  solid  surface  and  block  fluid  boundaries. 

The  time-marching  simulation  of  the  aeroelastic  responses  is  obtained  using  the  state  tran¬ 
sition  matrix  solution  from  t  to  t+At  of  the  state  variable  representation  of  the  decoupled 
modal  equations  [62,  65].  The  state  transition  matrix  based  scheme  is  optimal  in  the  sense 
that  it  is  derived  from  an  exact  solution  of  the  free  response  of  the  modal  equations.  The 
actual  scheme  uses  predictor/corrector  steps.  The  predictor  step  marches  the  structure  us¬ 
ing  the  solution  of  the  modal  equations  at  the  step  n  to  get  the  surface  deflection  at  the 
time  step  n+1.  This  provides  the  surface  shape  for  a  recomputation  of  the  fluid  mesh  and 
the  fluid  domain  solution  at  n+1.  After  a  solution  of  the  fluid  domain  involving  multiple 
subiterations,  the  corrector  step  then  solves  the  modal  equations  at  the  time  step  n+1  using 
the  averaged  generalized  forces  at  n  and  n+1. 

Because  the  CFD  and  CSM  meshes  usually  do  not  match  at  the  interface,  CFD/CSM  cou- 
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pling  requires  a  surface  spline  interpolation  between  the  two  domains.  The  interpolation 
of  CSM  mode  shapes  to  CFD  surface  grid  points  is  done  as  a  preprocessing  step.  Modal 
deflections  at  all  CFD  surface  grids  are  first  generated.  Modal  data  at  these  points  are  then 
segmented  based  on  the  splitting  of  the  flow  field  blocks.  Mode  shape  displacements  located 
at  CFD  surface  grid  points  of  each  segment  are  used  in  the  integration  of  the  generalized 
modal  forces  and  in  the  computation  of  the  deflected  surface.  The  final  surface  deformation 
at  each  time  step  is  a  linear  superposition  of  all  the  modal  deflections. 

The  following  is  an  account  of  our  theoretical  modeling  of  Schewe’s  experiment  on  transonic 
flutter  of  a  two-dimensional  supercritical  wing  with  an  NLR7301  airfoil  section  [70,  71,  75]. 
Figure  4.3  depicts  a  simplified  model  of  the  two-degree-of-freedom  test  set-up.  The  two- 
dimensional  wing  has  a  chord  length  of  0.3  m  (c  =  0.3  m)  and  a  span  of  1  m  (b  =  1  m). 
The  pitching  spring  and  the  plunging  spring  are  attached  to  the  same  c/4  position.  The 
corresponding  two-degree-of-freedom  equation  of  motion  of  the  set-up  reads 


'  mh  -sa  1  f  1  T  Dh  0  1  (  h\  \  Kh  0 

.  ~Sa  Ic/4  J  \  a  J  [  0  D*  J  l  «  J  +  [  0  K« 


where  m*  is  the  total  mass  (mh  =  26.64  kg),  Ic/4  is  the  mass  moment  of  inertia  about  c/4 
(/c/4  =  0.086  kg-m2),  sa  is  the  static  unbalance  (sQ  =  0.378  kg-m),  Dh  and  Da  are  the 
damping  factors  of  the  plunging  motion  (h)  and  the  pitching  motion  (a)  respectively  (Dh  = 
82.9  kg/s  and  Da  =  0.197  kg-m2/ (rad-s)),  Kh  and  Ka  are  the  stiffness  of  the  plunging  spring 
and  the  pitching  spring  respectively  (Kh  =  1.21  x  106  N/m  and  Ka  =  6.68  x  103  N-m/rad), 
and  L(t)  and  M(t)  are  the  aerodynamic  lift  and  moment  respectively  in  Newtons. 

The  aeroelastic  equations  and  the  CFD  grid  are  maintained  in  dimensional  form.  To  perform 
the  time-marching  CFD  computation  in  CFL3D  v6.0,  it  is  necessary  to  convert  equation  (1) 
into  modal  coordinates,  i.e., 


{  „  }  =  (4-2) 

where  q  is  the  modal  coordinate  and  (j)  is  the  modal  matrix  of  the  undamped  structure.  For 
this  numerical  example  we  have 


-0.1735  0.1004 
0.9277  3.403 


(4.3) 


Substituting  equation  (2)  into  equation  (1)  and  pre-multiplying  the  resulting  equation  by 
4>T  yields 
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where  Uh  and  uia  are  the  undamped  natural  frequencies  of  the  plunging  and  pitching  motions 
respectively  (u>h  —  205.4  rad/s  and  coa  =  299.3  rad/s),  and  (a  are  the  plunging  and  pitching 
damping  ratios  respectively  ((h  =  0.00649  and  (a  =  0.00521).  Note  that  the  off-diagonal 
terms  in  the  damping  matrix  are  assumed  to  be  zero  for  simplicity. 


4.3  Results  and  Discussion 

The  simulated  case  here  is  the  measurement  No.77  documented  in  [71].  As  mentioned  before, 
the  experimental  model  was  a  two-dimensional  supercritical  wing  with  NLR7301  section.  The 
chord  length  of  the  wing  was  0.3  m  and  the  angle  of  attack  was  1.28°.  The  experimental 
conditions  were  the  freestream  Mach  number  of  0.768  and  the  Reynolds  number  of  1.727  x 
106  based  on  the  chord  length.  A  transonic  LCO  in  two-degrees-of-freedom  was  found  at  the 
dynamic  pressure  of  0.126  bar.  The  corresponding  free-stream  velocity  was  254.7  m/s.  The 
total  pressure  was  0.45  bar. 

As  discussed  in  [60,  77],  because  of  the  relatively  large  chord  length  of  the  airfoil  with 
respect  to  the  wind  tunnel  test  section  (1  m  x  1  m ),  both  the  freestream  Mach  number  and 
the  angle  of  attack  need  to  be  corrected  to  take  into  account  wind  tunnel  wall  effects.  The 
criterion  used  in  [60,  77]  was  to  match  the  computed  to  the  measured  time-averaged  surface 
pressure  distribution. 


4.4  Time- Averaged  Surface  Pressure  Distribution 

An  Euler  computation  is  first  performed  on  a  C-type  grid  with  293  x  61  points.  The  best 
agreement  with  the  experimental  data  is  found  at  M  =  0.734  and  a  =  —0.25°.  Figure 
4.4  indicates  that  even  for  this  corrected  Mach  number  and  the  corrected  angle  of  attack, 
the  predicted  shock  strength  is  stronger  than  the  experimental  result  and  the  location  of 
the  shock  is  behind  the  measured  one.  It  seems  impossible  to  match  both  the  measured 
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x/c 


Figure  4.4:  Time- Averaged  Surface  Pressure  Distribution  (spring-off). 


Figure  4.5:  Time-Averaged  Surface  Pressure  Distribution  (spring-on). 


(a)  By  B-L  model  (b)  By  S-A  model 


Figure  4.6:  Entropy  Contours  (spring-on):  (a)  by  B-L-D-S  model;  (b)  by  S-A  model. 
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strength  and  location  of  the  shock  with  Euler  computation.  Two  viscous  computations  are 
then  performed  on  a  C-type  grid  with  293  x  93  points.  One  uses  the  Degani-Schiff  modified 
Baldwin-Lomax  turbulence  model  (B-L),  and  the  other  uses  Spalart-Allmaras  turbulence 
model  (S-A).  The  corrected  Mach  number  is  found  to  be  0.748  for  both  models  while  the 
corrected  angle  of  attack  is  -0.02°  for  the  B-L  model  and  0.15°  for  the  S-A  model.  It  is  found 
in  Figure  4.4  that  both  viscous  results  have  a  closer  agreement  with  the  experimental  data, 
especially  for  shock  strength  and  location,  clearly  showing  that  viscous  effects  are  important 
for  the  accurate  prediction  of  the  shock. 

The  above  computed  results  are  for  the  spring-off  condition.  With  taking  into  account 
the  contribution  from  the  structure  part,  the  corrected  angle  of  attack  becomes  -0.1°  for 
Euler  computation,  0.15°  for  the  B-L  model,  and  0.32°  for  the  S-A  model.  The  corrected 
Mach  number  does  not  change  for  both  Euler  computation  and  the  B-L  model  but  shifts 
to  0.75  for  the  S-A  model.  These  numbers  will  be  used  as  inputs  for  the  further  dynamic 
aeroelastic  computations. 

One  important  phenomenon  in  Figure  4.5  is  that  the  B-L  model  produces  a  pressure  or 
lift  loss  on  the  lower  surface  near  the  trailing  edge,  which  is  not  found  in  Figure  4.4,  the 
spring-off  condition.  In  order  to  understand  the  mechanism  behind  this  phenomenon,  Figure 
4.6  presents  the  entropy  contours  given  by  the  two  turbulence  models.  It  is  found  that  the 
B-L  model  seems  to  produce  spurious  vorticity  on  the  lower  surface  near  the  trailing  edge 
for  the  static  aeroelastic  computation. 


4.5  Performance  and  Time  Step  Convergence 

The  static  aeroelastic  case  is  used  to  compare  run  times  between  the  single  and  8  block/MPI 
parallel  computations.  Starting  from  a  steady  state  a  time-accurate  aeroelastic  solution  is 
marched  for  800  time  steps  at  5  multigrid  subiterations  each.  The  S-A  turbulence  model  is 
used  in  each  case. 

Table  1  gives  performance  data  for  the  various  grids,  run  modes,  and  code  configurations. 
The  first  two  rows  in  the  table  represent  sequential  computation  of  the  1  and  8  block  versions 
of  the  grid.  The  8-block  grid  takes  28%  less  CPU  time  than  the  1-block  grid,  apparently 
due  to  better  caching  of  the  smaller  sized  blocks.  The  parallel  computation  with  8  proces¬ 
sors  runs  7.8  times  faster  than  the  same  run  sequentially.  The  last  two  rows  present  the 
computational  effort  required  for  the  spring  analogy  smoothing  step  and  the  combined  TFI 
and  aeroelasticity.  Three  iterations  of  the  spring  analogy  scheme  add  about  1.5%  computing 
time  to  the  solution.  The  TFI  grid  movement  and  the  aeroelastic  computation  add  about 
10%  to  the  run  time. 

Table  2  shows  solution  behavior  for  the  single  block  grid  for  successive  time  step  sizes.  These 
are  computed  at  M0 0  =  0.753,  a  =  0.6°  and  a  dynamic  pressure  of  0.126  bar,  using  the  S-A 
model.  Initialization  is  accomplished  with  a  static  aeroelastic  solution,  followed  by  an  initial 
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Grid/Comp.  Type 

TFI 

Spr.  Anal.  Smoothing 

CPU  Time/Processor  (seconds) 

1  block/seq 

Yes 

Yes 

8800 

8  block/seq 

Yes 

Yes 

5430 

8  block/MPI 

Yes 

Yes 

690 

8  block/MPI 

Yes 

No 

680 

8  block/MPI* 

No 

No 

610 

*  also  without  aeroelas 

ticity 

Table  4.1:  Cost/Performance 


At 

0^  (mm) 

aa(Deg.) 

MHz) 

MHz) 

3.72 

34.5 

34.5 

0.040 

9.17 

3.24 

34.3 

34.3 

0.0125 

8.99 

3.17 

34.3 

34.3 

Table  4.2:  Time  Step  Convergence,  1  Block  Grid 


perturbation  of  the  dynamic  simulation  of  -0.00114  in  the  velocity  of  the  second  mode.  The 
time  step  sizes  give  80,  250,  and  800  time  steps  per  pitch/plunge  cycle.  Nine  subiterations 
per  time  step  are  used.  Columns  two  and  three  are  half  amplitudes  of  the  fully  developed 
LCO  plunge  and  pitch,  while  columns  four  and  five  are  the  plunge/pitch  frequencies.  Fre¬ 
quencies  and  amplitudes  are  computed  based  on  a  sampling  of  the  last  8-10  cycles  of  motion; 
from  this  sampling  the  data  appears  to  be  nearly  converged  at  the  smallest  time  step.  At 
the  largest  time  step,  even  after  100  cycles,  the  amplitude  slowly  continued  to  grow  while 
at  the  two  smaller  time  steps  the  amplitude  fully  converged  to  LCO.  In  all  of  the  remaining 
dynamic  computations,  the  smallest  time  step  size  with  nine  multigrid  subiterations  is  used. 

The  speed  increase  offered  by  computing  in  parallel  is  appealing.  There  are  trade  offs  of 
course  when  using  coarse  grain  parallelization.  Depending  on  the  block  splitting  and  prob¬ 
lem,  the  multiblock  computations  can  require  substantially  more  subiterations.  This  fact  is 
most  evident  in  the  problem  at  hand.  Figure  4.7  shows  a  comparison  between  the  same  LCO 
computed  on  a  single  flow  field  block  and  the  same  grid  divided  into  8  blocks.  The  8-block 
grid  is  shown  in  Figure  4.2.  By  viewing  the  right  hand  plots  in  that  figure  it  is  clear  that 
eventually  the  amplitude  and  frequency  of  the  multiblock  pitch  and  plunge  had  settled  out 
and  converged  to  values  virtually  identical  to  that  of  the  single  block  grid  (a^  =  8.87  mm, 
aa  —  3.13°  ,  Uh  =  ua  =  34.3  Hz).  This  would  not  be  the  case  if  the  multiblock  aeroelastic 
coupling  and  integration  were  not  consistent  with  that  of  the  single  block  configuration.  Yet 
as  shown  in  Figure  4.7  the  time  at  which  the  LCO  can  be  considered  converged  is  different 
between  the  two  computations.  This  accounts  for  the  different  values  of  to  for  the  1  block 
and  8  block  computations.  As  seen  in  the  left  hand  plots  in  Figure  4.7,  the  LCO  computed 
with  8  flow  field  blocks  takes  much  longer  to  develop  than  does  an  identical  LCO  computed 
on  a  single  flow  field  block. 
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(a)  Plunging  motion 
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(b)  Pitching  motion 


Figure  4.7:  Comparison  of  Single/Multiblock  Grid  Computations  (N-S,  =  0.753,  a  = 

0.6°)  t0=4.5s  (8  blocks):  (a)  Plunge  Motion;  (b)  Pitch  Motion. 


4.6  Effect  of  Viscosity  and  Turbulence  Modeling 

Figure  4.8  compares  LCO  predicted  by  Euler  and  Navier-Stokes  methods.  It  is  found  that 
although  both  inviscid  and  viscous  computations  are  able  to  predict  LCO,  the  amplitude 
of  LCO  given  by  viscous  computations  is  less  than  1/3  of  that  predicted  by  inviscid  com¬ 
putation,  much  closer  to  the  experimental  data  (plunging  and  pitching  amplitudes  are  0.65 
mm  and  0.18°  respectively).  Therefore,  it  is  important  to  include  the  viscous  terms  in  LCO 
computations,  which  effect  is  to  limit  the  amplitude  of  the  flutter  LCO. 

Also  shown  in  Figure  4.8  is  the  difference  between  LCO  predicted  by  the  two  turbulence 
models.  It  is  found  that  the  B-L  model  reaches  large  amplitude  limit  cycle  quite  rapidly. 
Plunge  and  pitch  amplitudes  of  the  B-L  result  are  around  4  mm  and  2°  respectively.  The 
predicted  frequency  is  32.2  Hz.  On  the  other  hand,  the  S-A  model  reaches  LCO  much  later. 
Its  plunge  and  pitch  amplitudes  are  larger  than  those  given  by  the  B-L  model,  around  9  mm 
and  3°  respectively.  Furthermore,  the  B-L  LCO  is  a  random/chaotic  output  while  the  S-A 
LCO  is  a  prefect  sine- wave  output.  It  is  clear  that  the  turbulence  model  significantly  alters 
the  nature  of  the  solution. 

In  order  to  understand  the  mechanism  behind  the  differences  caused  by  the  two  turbulence 
models,  Figure  4.9  further  presents  the  predicted  density  contours  at  8  time  steps  uniformly 
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(a)  By  B-L  model 


Figure  4.9:  Predicted  Density  Contours: 


(b)  By  S-A  model 
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distributed  in  the  last  cycle.  The  B-L  results  presented  at  the  left  hand  side  clearly  show 
multishocks  with  large  execusions  possibly  caused  by  spurious  vorticity  separation.  This  is 
why  the  B-L  LCO  exhibits  chaotic  behavior.  On  the  other  hand,  no  visible  vorticity  sep¬ 
aration  is  found  in  the  S-A  results  presented  at  the  right  hand  side.  As  a  result,  the  flow 
pattern  is  much  more  regular. 


4.7  Effect  of  Perturbation  Size 

The  above  LCO  computations  are  performed  with  a  small  initial  perturbation  (-0.0012  in 
the  velocity  of  the  second  mode).  In  order  to  investigate  the  effect  of  initial  perturbation 
size,  the  earlier  simulations  are  repeated  with  the  two  larger  initial  perturbations.  The  one 
referred  to  as  the  medium  perturbation  in  Figures  4.10  and  4.11  has  5  times  of  the  original 
perturbation  size  (-0.006  in  the  velocity  of  the  second  mode),  and  that  referred  to  as  the 
large  perturbation  has  10  times  of  the  original  perturbation  size  (-0.012  in  the  velocity  of 
the  second  mode). 

Results  given  by  the  B-L  model  are  shown  in  Figure  4.10.  It  is  found  that  although  the  LCO 
amplitudes  are  in  the  same  level  for  small  and  medium  perturbations,  the  LCO  amplitudes  for 
large  perturbation  is  significantly  reduced.  This  is  contradicted  to  common  senses.  However, 
this  phenomenon  can  even  be  more  clearly  seen  in  Figure  4.11,  which  are  S-A  results.  From 
small  perturbation  to  medium  perturbation,  LCO  occurs  earlier  and  the  amplitudes  of  both 
plunging  and  pitching  modes  are  slightly  reduced.  Further  increasing  the  perturbation  size 
by  two,  a  converged  solution  is  achieved.  Larger  perturbation  sizes  correspond  to  stronger 
nonlinear  aerodynamic  effects.  So,  the  above  results  mean  that  nonlinear  aerodynamic  effects 
tend  to  stabilize  the  flutter. 


4.8  Conclusions 

Stronger  nonlinear  aerodynamic  effects,  including  adding  the  viscous  terms  and  increasing 
perturbation  size,  tend  to  reduce  LCO  amplitudes  and  even  stabilize  the  dynamic  solution. 
Viscous  LCO  computation  is  turbulence  model  dependent.  Specifically,  B-L  model  with 
Degani-Schiff  modification  results  in  multishocks  with  large  execusions  possibly  caused  by 
spurious  vorticity  generation. 
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(b)  Pitching  mode 


Figure  4.10:  Effect  of  Initial  Perturbation  Sizes:  (a)  Plunge  Motion;  (b)  Pitch  Motion 


(b)  Pitching  mode 


Figure  4.11:  Density  Contours  Showing  Effect  of  Initial  Perturbation  Sizes:  (a)  1st  Cycle; 
(b)  2nd  Cycle;  (c)  3rd  Cycle;  (d)  Final  LCO  Cycle. 


Chapter  5 


Theoretical  Background  of  the 
Harmonic  Balance  Method  in 
Three-Dimensions 

Summary 

In  order  to  facilitate  the  three-dimensional  work  in  the  following  chapters.  Theoretical  back¬ 
ground  of  the  Harmonic  Balance  Method  in  three  dimensions  will  be  described  in  this  chapter. 
The  sections  that  follow  include: 

•  Three  Dimensional  Navier-Stokes  and  Euler  Equations 

•  Harmonic  Balance 

•  Time-Linearized  Flow/Three  Dimensional  Flow  Solver 

•  Reduced  Order  Modeling  of  Time-Linearized  Aerodynamic  Models 

•  Reduced  Order  Aeroelastic  Model  for  Multi- Degree- Of- Freedom  Structures 

•  Structural/ Compatible  Reduced  Order  Aeroelastic  Model 

5.1  Three-Dimensional  Navier-Stokes  and  Euler  Equa¬ 
tions 

We  start  with  the  top-down  approach: 

The  3-D  Reynolds-Averaged  Navier-Stokes  Equations  read: 
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/  \ 

p 

f  \ 

pu 

f  \ 

pv 

f  \ 

pw 

pu 

pu 2  +p 

puv 

puw 

pv 

>  ,F-  < 

puv 

>,G=  < 

pv2  +p 

>,H  =  < 

pvw 

>, 

pw 

puw 

pvw 

pw2  +  p 

,  Pe  , 

<  puh  j 

,  Pvh  , 

<  pwh  J 

'  0 

r  0 

r  0  ' 

Txx 

Tyx 

7~ZX 

Txy 

II 

Tyy 

II 

TZy  > 

T~xz 

ryz 

T~zz 

w  Txx'U>  “1“  TXyV  “1“  Txz'W  J 

<  TyxU  +  TyyV  +  TyZW 

,  TZXU  +  TzyV  +  TZZW  j 

and  p,  p,  (u,  v,  w),  e,  h  are  the  flow  density,  pressure,  velocity,  energy  and  enthalpy,  respec¬ 
tively.  Txx,rxy,Txz,  •,  etc  are  the  Reynolds  stress  tensors. 

In  what  follows,  we  will  reduce  Eq  (5.1)  to  3-D  and  2-D  Euler  equations  for  the  conve¬ 
nience  of  elucidating  the  frequency-domain  POD/ROM  approach. 


5.2  Harmonic  Balance 

In  Phase  I,  the  importance  of  aerodynamic  nonlinearities  on  LCO  was  assessed  using  a  novel 
two-dimensional  harmonic  balance  (H.B.)  technique.  With  this  approach,  the  unsteady  flow 
variables  can  be  represented  by  a  Fourier  Series  in  time  with  spatially  varying  coefficients. 
This  assumption  leads  to  a  set  of  harmonic  balance  Euler  equations,  which  can  be  solved  ef¬ 
ficiently  using  conventional  CFD  methods  including  time  marching  with  local  time  stepping 
and  multi-grid  acceleration.  The  two-dimensional  Euler  equations  are: 


dU  dF  dG 
dt  dx  partialy 


(5.2) 


where  the  flux  vectors  U,  F  and  G  are  given  by 
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Next,  representing  conservation  prime  variables  by  sum  of  harmonics  of  fundamental  fre¬ 
quency  gives: 
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(5.4) 


p(x,y,t)  =  ^  RJx,  y)e'unt,pa(x,  y,  t)  =  UJx,  y)eM 

n  n 

pv(x,y,t)  =  ^2vn(x,y)eiunt,pe(x,y,t)  =  ^En{x,y)eiwnt 

n  n 


Requiring  each  frequency  component  to  vanish  independently  (harmonic  balance)  and  col¬ 
lecting  terms  of  like  harmonics  results  in 


where 


and 


duty)  dF(V)  dGiy) 

dt  dx  dy 


=  0 


V  —  {. . .  Rq,  Uq,  Vq,  Eo,  R+ i,  U+i,  V+i,  E+ 1, . .  .}T 


(5.5) 


dJJ 

—  iv{'  •  ’0  •  Rq,  0  •  Uo,  0  •  Vo,  0  •  Eo,  +1  •  R+x,  +1  •  U+ i,  +1  •  V+i,  +1  •  E+  i}T  (5.6) 

By  adding  a  ’pseudo- time’  term,  Eq  (5.5)  can  be  solved  by  a  conventional  CFD  solution,  i.e. , 

ay  dF(v)  dG(v)  eu(v) 

d7  +  ~a^T  +  ~dT  +  ~sr  =  0  (5-7) 


Note  that  harmonic  balance  equations  are  similar  in  form  to  original  Euler  equations.  Thus, 
existing  Euler  codes  can  be  modified  to  solve  H.B.  equations.  A  two-step  Lax-Wendroff 
scheme  will  be  used.  Also,  since  only  ’’steady-state”  solution  is  desired,  one  can  use  local 
time  stepping,  multiple-grid  acceleration  techniques  and  residual  smoothing  to  speed  conver¬ 
gence.  The  method  is  computationally  very  efficient,  at  least  one  to  two  orders-of-magnitude 
faster  than  nonlinear  time-domain  CFD  simulations. 


In  Phase  II,  we  will  use  the  two-dimensional,  nonlinear,  harmonic  balance  aerodynamic 
code  and  add  viscosity  terms,  i.e.,  we  will  develop  a  Navier-Stokes  version  of  this  analysis 
suitable  for  the  analysis  of  oscillating  airfoils.  In  addition,  we  will  use  this  code  to  perform 
LCO  analysis  and  prediction  for  simple  structural  models,  e.g.  plunge  and  pitch  of  a  typical 
section  for  conventional  and  supercritical  airfoils.  In  particular,  we  will  perform  correlation 
studies  with  the  DLR  experimental  data  of  the  NLR  7301  LCO  case. 

In  Phase  I,  a  three-dimensional  time-linearized  unsteady  aerodynamic  analysis  was  devel¬ 
oped.  In  Phase  II,  this  code  will  be  converted  to  a  harmonic  balance  analysis  of  three- 
dimensional  inviscid  transonic  flows.  We  will  use  this  three-dimensional  analysis  to  perform 
LCO  analyses  of  simple  linear  wing  structural  models.  We  will  perform  LCO  correlation 
studies  with  a  business  jet  wing  or  other  wing  to  be  determined. 
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In  passing,  we  note  the  following  when  using  aeroelastic  eigenmodes  for  LCO  prediction. 
The  harmonic  balance  method  has  proven  to  be  a  very  effective  method  for  determining 
limit  cycle  oscillations  (LCO)  and  other  nonlinear  responses  of  aeroelastic  systems.  It  is 
often  more  computationally  efficient  and  gives  greater  physical  insight  than,  for  example, 
time  marching  simulations.  On  the  other  hand  it  has  some  limitations,  e.g.  it  does  not  allow 
an  investigation  of  the  possible  dependence  of  the  response  on  the  initial  conditions.  Nor 
does  it  allow  a  precise  assessment  of  non-periodic  and  chaotic  motions.  At  best  the  harmonic 
balance  method  can  only  approximate  a  non-periodic  motion.  In  practice  these  limitations 
are  not  often  serious  drawbacks;  however  these  observations  do  suggest  that  future  work 
may  well  involve  a  balanced  consideration  of  both  time  simulations  and  frequency  domain 
(harmonic  balance)  studies. 

Another  consideration  is  the  number  of  harmonics  to  be  retained  in  a  harmonic  balance 
analysis.  If  there  is  a  dominant  harmonic,  a  single  harmonic  may  suffice,  although  it  will  be 
important  to  assess  the  effects  of  higher  harmonics  on  the  fundamental  harmonic  per  se.  And 
finally  the  algebraic  complexity  of  using  the  method  for  multi-modal  systems  may  quickly 
get  out  of  hand  unless  special  measures  are  taken  to  deal  with  this  issue.  We  will  return 
to  this  latter  issue  as,  even  with  reduced  order  aerodynamic  modeling  of  the  structure  and 
fluid,  the  order  of  the  reduced  model  will  be  larger  than  that  often  encountered  in  classical 
harmonic  balance  analyses,  i.e.,  on  the  order  of  several  tens  of  states. 


5.3  Time-Linearized  Flow  Three-Dimensional  Flow  Solver 

In  Phase  I,  we  developed  a  three-dimensional  time-linearized  Euler  analysis  of  unsteady  flow 
about  a  wing.  We  start  with  the  assumption  that  the  unsteady  flow  u(x,  y,  z,  t )  may  be 
expanded  in  a  perturbation  series  of  the  form 


u(x,  y,  z,t )  =  U(®,  y,  z)eiwt 


(5.8) 


where  U (x,ytz,t)  represents  the  steady  background  flow,  and  u (x,y,z,t)  is  the  complex 
amplitude  of  the  unsteady  small-disturbance  flow  arising  from  the  wing  vibration,  which 
vibrates  at  frequency  u.  Substituting  Eq  (5.8)  into  the  nonlinear  three-dimensional  Euler 
equations,  and  expanding  the  result  in  a  perturbation  series  in  the  small-disturbance  quan¬ 
tities,  one  finds  that  to  zeroth  and  first  order  the  governing  equations  reads  respectively 


<9F(U)  dG(U)  dH(U) 
dx  dy  dz 


d 

iuu+  — 
ox 


(5.9) 

(5.10) 


In  a  conventional  time-linearized  analysis,  Eq  (5.9)  and  Eq  (5.10)  are  solved  using  pseudo 
time  marching  and  standard  CFD  techniques.  Note  that  since  time  does  not  appear  explic¬ 
itly  in  either  equation,  they  both  may  be  solved  using  computational  acceleration  techniques 
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such  as  multiple  grid  acceleration,  local  time  stepping,  and  residual  smoothing.  The  result  is 
that  these  two  equations  can  be  solved  at  a  fraction  of  the  cost  of  solving  the  full  nonlinear 
Euler  equations  using  conventional  time  marching. 

In  Phase  I,  a  simple  clean  wing  model  was  developed.  This  CFD  analysis  used  a  simple 
grid  structure  that  limited  the  approach  to  clean  wings.  In  Phase  II,  this  model  will  be  ex¬ 
tended  to  allow  more  complicated  grid  topologies,  which  will  allow  us  to  analyze  multi-body 
configurations,  e.g.  wing/store. 

5.4  Reduced  Order  Modeling  of  Time-Linearized  Aero¬ 
dynamic  Models 

Having  developed  a  three-dimensional  time-linearized  flow  solver,  we  next  consider  the  re¬ 
duction  of  that  model  using  proper  orthogonal  decomposition  techniques  (POD).  The  idea 
behind  the  frequency  domain  proper  orthogonal  decomposition  is  a  simple  one.  We  first 
calculate  the  small-disturbance  response  of  the  aerodynamic  system  at  M  different  combina¬ 
tions  of  frequency  and  excitation.  The  solutions,  also  called  ’’snapshots,”  are  denoted  by  qm 
for  m=l,2„M.  These  snapshots  are  then  linearly  combined  to  for  a  smaller  number  of  basis 
vectors  <f>k  for  k=l,2„K  where  KjM.  In  other  words, 

4>k  =  Q  vkforJc  =  1, 2, ...,  K  (5.11) 

where  the  rath  column  of  Q  is  just  qm.  In  the  proper  orthogonal  decomposition  technique, 
the  vectors  v*  are  found  by  solving  a  small  eigenvalue  problem  of  the  form 

Q^Qvfc  =  (5.12) 


Only  the  eigenvectors  vfc  with  the  largest  eigenvalues  Xk  are  used  to  form  basis  vectors  de¬ 
fined  by  Eq  (5.11).  QH  is  the  Hermitian  of  Q. 

Having  computed  the  POD  basis  vectors,  we  assume  that  they  will  provide  a  useful  ba¬ 
sis  for  computing  the  unsteady  solution  at  some  other  frequency  and/or  external  excitation 
than  was  use  to  general  the  original  snapshots.  Thus,  we  let 

q  =  (5.13) 

where  $  is  an  NxK  matrix  whose  kth  column  is  simply  the  basis  vector  4>k,  and  £  is  a  vector 
of  aerodynamics  state  variables. 

We  note  that  when  discretized,  Eq  (5.10)  has  the  form 

Aq  =  A0q  +  icuAiq  =  b0  +  iwbi  (5. 14) 
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where  Ao  and  Ai  are  independent  of  the  excitation  frequency  co,  and  are  purely  real,  and  q 
is  the  vector  containing  u  at  all  the  nodes  of  the  CFD  grid. 

Substitution  of  Eq  (5.11)  into  Eq  (5.14),  and  projection  of  the  error  onto  the  space  spanned 
by  the  basis  vectors  gives 

=  A  r£  =  (5.15) 


Finally,  the  reduced-order  aerodynamic  matrix  AR,  which  is  much  smaller  than  the  original 
aerodynamic  matrix  A,  is  factored  using  LU  decomposition,  and  Eq  (5.15)  is  solved  for  the 
unknown  aerodynamic  state  variables.  This  step  is  computationally  very  efficient  because 
the  reduced-order  aerodynamic  matrix  is  small,  sometimes  as  small  as  10x10,  but  rarely 
larger  than  100x100.  The  major  expense  in  constructing  the  reduced-order  aerodynamic 
model  is  the  computation  of  the  snapshots;  the  cost  of  finding  the  basis  vectors  and  solution 
to  Eq  (5.15)  is  negligible  by  comparison. 

Having  described  the  basic  reduced-order  modeling  technique,  we  next  describe  how  to  incor¬ 
porate  an  aerodynamic  reduced-order  model  into  an  aeroelastic  model  of  flutter.  Consider, 
for  example,  a  two  degree-of-freedom  structural  dynamic  model  of  a  typical  section.  The 
governing  equations  of  motion  are  of  the  form 

Mh  +  Kh  =  fwhereh.  =  { h ,  a}T  (5.16) 

and  h  and  a  are  the  plunging  and  pitching  degrees  of  freedom  of  the  typical  section.  M  and 
K  are  the  mass  and  stiffness  matrices,  respectively. 

Note  that  the  aerodynamic  force  vector  f  is  obtained  from  integrals  involving  the  pres¬ 
sure  at  the  surface  of  the  airfoil.  When  discretized,  these  integrals  may  be  expressed  as 

f  =  Cq  (5.17) 


where  C  is  a  sparse  2 xN  matrix.  Similarly,  for  the  case  of  airfoil  vibration,  the  vector  b  on 
the  right-hand  side  of  Eq  (5.14)  can  be  expressed  as 

b  =  bo  +  icubi  =  Boh  +  icuBih  (5.18) 

where  now  we  have  made  the  assumption  that  the  airfoil  motion  is  harmonic  in  time  al¬ 
though  u  may  be  complex).  For  large  CFD  models,  finding  the  eigenvalues  is  prohibitively 
expensive.  To  reduce  the  size  of  the  model,  we  again  assume  that  the  number  of  aerodynamic 
states  can  be  reduced  using  Eq  (5.13).  Again,  projecting  the  error  of  the  aerodynamic  equa¬ 
tions  onto  the  space  spanned  by  the  aeroelastic  basis  vectors  gives  the  desired  reduced-order 
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aeroelastic  model,  i.e. 


■  -$hb„  -$hb!  i  r  )  aRi  o  o  i  r  r 

0  0  I  h  Uiw  0  -I  0  h=0  (5.19) 

-C$  K  0  h  0  0  M  h 

L  J  V  L  J  L  J 

5.5  Reduced-Order  Aeroelastic  Model  for  Multi-Degree- 
Of-Freedom  Structures 

In  Eq  (5.19),  we  have  derived  a  reduced-order  aeroelastic  model  for  a  two  degree  of  freedom 
(2  d.o.f.)  structure,  i.e.,  h  =  h,aT.  For  multi-degree-of-freedom  structures,  h  can  also  be 
interpreted  as  the  displacements  of  each  degree  of  freedom  in  the  structural  model.  In  this 
situation,  the  size  of  Eq  (5.19)  becomes  K  +  2 n,  where  n  is  the  number  of  structural  degrees 
of  freedom.  For  a  realistic  wing  structural  model,  n  can  be  on  the  order  of  thousands, 
rendering  Eq  (5.19)  a  very  large  size  eigenvalue  problem.  Solving  such  a  large  eigenvalue 
problem  would  practically  be  impossible. 

5.5.1  The  Modal  Approach 

To  circumvent  this  problem,  we  introduce  the  so-called  ”modalapproach"  to  Eq  (5.19).  The 
modal  approach  approximates  h  by  the  superposition  of  the  low-order  structural  modes,  i.e., 

h  =  ’F ar]  (5.20) 

where  \I/a  is  the  modal  matrix  whose  columns  contain  the  modal  data  of  the  low-order  struc¬ 
tural  modes  and  77  are  is  the  generalized  coordinate  vector.  Since  the  magnitude  of  the 
modes  can  be  arbitrary,  they  are  usually  normalized  by  the  square  root  of  their  respective 
generalized  masses,  giving  a  unit  generalized  mass  matrix.  The  justification  for  using  the 
modal  approach  is  based  on  the  premise  that  most  of  the  aeroelastic  responses  are  dominated 
by  the  lower-order  structural  modes.  Usually,  for  a  single  wing  structure,  the  lowest  ten  (10) 
structural  modes  axe  sufficient  to  accurately  represent  h.  Substituting  Eq  (5.20)  into  Eq 
(5.19)  and  pre-multiplying  the  second  and  third  rows  of  Eq  (5.19)  by  yields 

$HB01Fa  -$HB!$a  1  (  n  r  ARl  0  0  1  £ 

0  [A]  (  >+*»  0  [a]  0  <  V  =  0(5.21) 

K]  [2w„,c]  J  i  J  0  0  I J  l  4  J 

where  uini  and  m*  are  the  natural  frequency  and  the  generalized  mass  of  the  ith  structural 
mode,  respectively.  Now,  the  size  of  Eq  (5.21)  becomes  K  +  2m,  where  m  is  the  number  of 
structural  modes.  Because  m  is  much  less  than  the  number  of  structural  d.o.f.  n,  the  size 
of  the  aeroelastic  system  is  substantially  reduced. 


Aro  — 
0 
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5.5.2  The  Structural/Modal  Damping 

Note  that  a  structural/modal  damping  matrix  diag  [. . .  2u;ni£ . . .]  has  been  added  to  Eq  (5.21). 
With  this  added  matrix,  Eq  (5.21)  will  facilitate  our  subsequent  study  of  LCO.  In  an  earlier 
LCO  study,  ZONA  has  suggested  that  the  structural  damping  could  play  a  decisive  role 
in  LCO  for  a  wing/store  system.  For  LCO  study  using  the  proposed  method,  one  could 
alternatively  use  Eq  (5.21)  to  delineate  the  effects  with  and  without  structural  damping. 
Also  note  that  the  modal  approach  in  fact  increases  the  sparseness  of  the  matrices  in  Eq 
(5.21).  Thus,  the  eigenvalue  problem  of  Eq  (5.21)  can  be  solved  more  efficiently  by  a  sparse 
eigensolver. 

5.6  Structural  Compatible  Reduced-Order  Aeroelastic 
Model 

One  of  the  technical  issues  involved  in  the  CFD  aeroelastic  computations  is  the  displacement 
transferal  from  the  structural  finite  element  grid  to  the  CFD  surface  grid.  This  technical 
issue  arises  from  the  problem  where  the  CFD  surface  grid  and  the  structural  finite  element 
(FEM)  grid  are  considerably  different;  in  their  locations  and  densities.  Solving  such  a  dis¬ 
placement  transferal  problem  of  complex  configuration  such  as  whole  aircraft  with  external 
stores  is  by  no  means  a  trivial  task. 

In  March  1999,  ZONA  received  a  two-year  contract  from  NASA/Langley  to  develop  a  Bound¬ 
ary  Element  Method  (BEM),  called  the  BEM  Solver,  for  the  data  transferal  between  the  FEM 
grid  and  the  CFD  surface  grid  (see  Section  6.1,  ZONA’s  Related  Work).  By  formulating  the 
data  transferal  problem  as  an  equivalent  solid  mechanics  problem,  the  BEM  Solver  generates 
a  universal  spline  matrix  [S]  which  relates  the  displacements  at  the  FEM  grid  to  the  CFD 
grid  such  that 

Va  =  S*s  (5.22) 

where  and  are,  respectively,  the  modal  matrix  at  the  CFD  grid  and  at  the  FEM 
grid.  With  the  universal  spline  matrix  [S]  at  hand,  converting  Eq  (5.21)  to  a  structural- 
compatible  reduced-order  aeroelastic  model  is  straightforward.  Substituting  Eq  (5.22)  into 
Eq  (5.21)  gives: 

Ar,,  -$hB0S*s  -$hBxS*s  1  f  \  ARl  0  0 1  r  £  \ 

0  0  [i]  \  v\+™  0  [i]  0  <  V  >=(0-23) 

_  -*jsTc<r>  [<]  [2u>niC]  J  U  J  o  o  i  J  U 

Finally,  we  arrive  at  the  structural-compatible,  modal-based,  reduced-order  aeroelastic  model, 
Eq  (5.23).  The  size  of  this  model  is  K  +  2m,  where  m  is  the  number  of  low-order  structural 
modes  and  is  usually  on  the  order  of  ten  for  a  single  wing  structure.  Eq  (5.23)  contains  two 
very  sparse  matrices  that  can  be  solved  efficiently  by  a  sparse  eigenvalue  solver. 
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5.7  Modeling  LCO  of  High  Degree-of-Freedom  Nonlin¬ 
ear  Systems 

The  nonlinear  system  that  we  wish  to  model  can  for  the  most  part  be  modeled  as  quasi-linear 
nonlinear  vector  equations  of  the  form 

<wMq  +  Nq  =  0  (5.24) 


where  q  is  a  very  large  matrix  arising  containing  the  Fourier  coefficients  of  the  unknowns  in  a 
harmonic  balance  analysis,  M  is  the  a  constant  matrix,  and  N  is  a  nonlinear  vector  operator 
arising  from  the  harmonic  balance  analysis.  Here  q  would  contains  the  unknown  flow  solution 
in  the  three-dimensional  harmonic  balance  of  the  three-dimensional  flow  field,  and  also  the 
structural  dynamic  modes.  Thus,  the  system  may  contain  hundreds  of  thousands  of  degrees 
of  freedom.  This  system  of  equations  is  solved  using  pseudo-time  time  marching.  Thus,  Eq 

(5.24)  is  solved  by  marching  the  equation 

^-Mu;  Mq  +  Nq  =  0  (5.25) 

in  time  until  a  steady  state  is  reached.  However,  if  u)  is  not  known  accurately,  then  Eq 

(5.25)  will  not  converge,  but  will  itself  go  into  a  mathematical  limit  cycle.  One  can  show, 
however,  that  the  solution  q  will  be  nearly  correct.  The  solution  can  be  improved  by  first 
computed  a  better  estimate  for  u  using  a  nonlinear  Rayleigh  quotient,  i.e. 


O’  « 


qHN(q) 

qHMHM(q) 


(5.26) 


Eq  (5.25)  can  then  be  marched  again  to  improve  the  estimate  of  q.  This  process  can  be 
repeated  until  convergence.  The  result  is  a  description  of  the  LCO  behavior  of  a  very  high 
dimensional  system,  i.e.  a  nonlinear  CFD  model  coupled  to  a  linear  or  nonlinear  structural 
model.  If  Phase  II,  this  technique  will  be  applied  to  the  harmonic  balance  flow  solver  coupled 
to  a  linear  and/or  nonlinear  structure. 
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Chapter  6 

Development  of  the  Structural/FEM 
Models  and  CSD/CFD  Grid 
Interfacing  for  Candidate  Wings  for 
LCO  Investigation 


Summary 

•  In  support  of  the  CFD  activities  for  3D  LCO  investigation,  several  wing  planforms 
were  selected  as  candidates  for  LCO  case  studies.  The  selected  wing  and  wing /store 
candidates  were: 


1.  the  MAVRIC-I  Business  Jet  Wing  (Fig.  6.1) 

2.  the  ARW-2  Wing  (Fig.  6.2) 

3.  the  AGARD  Standard  445-6  Wing  (Fig.  6.3),  and 

4-  the  F-16  Wing  and  Wing/Store  System.  (Fig.  6.4) 

All  wings  under  consideration  are  known  to  have  occurrance  of  LCO  in  the  lower  tran¬ 
sonic  regime  by  Wind-Tunnel  testing  (at  TDT  of  NAS A-Langley). 

•  The  Configurations  of  these  wing  planforms  are  shown  below. 

•  One  supporting  efforts  to  the  follow-on  3D  LCO  investigations  through  out  the  Phase 
II  period  is  to  supply  structural/FEM  models  and  to  develop  CSD/CFD  interfacing 
models.  For  the  selected  wing-planform  candidates.  In  what  follows,  we  will  show 
mode  shapes  due  to  structural/FEM  versus  deformed  CFD  meshes  for  various  wing 
candidates  as  a  result  of  these  efforts. 

•  The  finally  selected  wing-planform  for  3D  LCO  case  studies  turned  out  to  the  AGARD 
Standard  445-6  wing  and  the  F-16  wing  and  wing/store  system.  Due  to  the  accessibility 
and  availability  of  their  test  data,  the  MAVRIC-I  wing  and  ARW-2  wing  were  not 
selected  for  the  preliminary  test  cases  for  Phase  II.  The  LCO  results  of  the  AGARD 
445.6  and  the  F-16  wing/store  system  will  be  presented  in  the  next  two  chapters. 
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Pressure  Transducer  Rows 
(28  per  row,  18  upper,  10  lower) 


Figure  6.1:  Instrumentation  Layout  for  Refurbished  MAVRIC-I  Business  Jet  Wing  Model 


ZY 
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Figure  6.3:  AGARD  Standard  445.6  wing,  Showing  Dimensions:  AR  =  4.0,  A=45  deg,  A=0.6, 
NACA  65A004  airfoil  section  (taken  from  AGARD  Report  765  by  Dr.  Carson  Yates.) 


Figure  6.4:  Wing  Aerodynamic  and  Structural  Model 


Mode 

Natural  Frequencies  (Hz) 

Generalized  Mass  (lbf-in/s2) 

Generalized  Stiffness  (lbf-in) 

1 

4.083394E+00 

1.000000E+00 

6.582675E+02 

2 

1.397037E+01 

1.000000E+00 

7.705050E+03 

3 

3.154093E+01 

1.000000E+00 

3.927433E+04 

4 

3.199304E+01 

1.000000E+00 

4.040833E+04 

5 

5.811350E+01 

1.000000E+00 

1.333257E+05 

6 

5.879132E+01 

1.000000E+00 

1.364540E+05 

7 

8.823108E+01 

1.000000E+00 

3.073286E+05 

8 

9.221439E+01 

l.OOOOOOE+OO 

3.357045E+05 

Table  6.1:  Natural  Frequencies,  Generalized  Masses,  and  Generalized  Stiffness  of  the 
MAVRIC-I  Business  Jet  Wing. 


6.1  Deforming  CFD  Meshes  of  the  MAVRIC-I  Busi¬ 
ness  Jet  Wing 

The  MAVRIC-I  business  jet  wing  is  an  ongoing  wind  tunnel  test  program  of  the  NASA/Langley 
aeroelasticity  branch.  The  model  is  constructed  with  a  simple  stepped  aluminum  plate  pro¬ 
viding  the  wing  stiffness  and  fitted  with  end-grain  balsa  wood  to  provide  the  wing  contour. 
Previous  experiment  has  shown  that  in  the  0.8  to  0.9  Mach  number  range,  the  model  motions 
were  predominantly  in  the  first  wing-bending  mode  and  exhibited  Limit  Cycle  Oscillation 
(LCO). 

The  model  has  been  retested  in  the  Transonic  Dynamic  Tunnel  (TDT)  as  the  MAVRIC- 
I.  Figure  6.1  indicates  the  location  of  instrumentation  that  has  been  added  to  the  model. 
Three  chords  of  unsteady  pressure  transducers  are  installed  at  span  stations  0.22,  0.63  and 
0.87.  Each  chord  has  28  upper  surface  and  18  lower  surface  close-mounted  transducers. 
Eight  accelerometers  are  mounted  along  the  leading  and  trailing  edge  of  the  wing  and  bend¬ 
ing/torsion  strain  gages  are  installed  at  the  root.  The  intent  of  the  retest  is  to  obtain 
unsteady  pressure  and  wing  response  data  under  conditions  of  LCO  in  order  to  validate 
CFD  codes  for  such  conditions. 

Figure  6.5  shows  a  structural  finite  element  model  of  the  MAVRIC-I  Business  Jet  Wing  that 
was  provided  by  the  NASA/Langley  aeroelasticity  branch.  This  structural  finite  element 
model  consists  of  613  grid  points  and  551  plate  elements  whose  first  8  natural  frequencies, 
generalized  masses,  and  generalized  stiffness  are  presented  in  Table  1. 

Figure  6.6  presents  a  surface  CFD  mesh  of  the  MAVRIC-I  business  wing  that  consists  of 
33x121  grid  points.  To  generate  the  deformed  CFD  surface  meshes  of  the  8  structural  modes, 
we  employed  the  ”BEM  Solver”  which  was  developed  by  ZONA  under  a  two-year  contrac¬ 
tual  support  of  NASA/Langley.  The  BEM  solver  utilizes  the  structural  boundary  element 
method  to  generate  a  universal  spline  matrix  that  can  transfer  the  displacement  computed 
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Figure  6.6:  CFD  Surface  Mesh  of  the  MAVRIC-I  Business  Jet  Wing. 


Mode 

Natural  Frequencies  (Hz) 

Generalized  Mass  (lbf-in/s2) 

Generalized  Stiffness  (lbf-in) 

1 

9.943 

1.0 

2.49E+03 

2 

29.85 

1.0 

3.159E+04 

3 

32.06 

1.0 

4.059E+04 

4 

61.13 

1.0 

1.475E+05 

5 

63.82 

1.0 

1.6079E+05 

6 

99.16 

1.0 

3.88E+05 

7 

105.718 

1.0 

4.412E+05 

8 

115.32 

1.0 

5.251E+05  ' 

Table  6.2:  Natural  Frequencies,  Generalized  Masses,  and  Generalized  Stiffness  of  the  ARW-2 
Wing. 


at  the  structural  grid  to  the  CFD  grid.  These  deformed  structural  FEM  models  computed 
by  MSC/NASTRAN  and  8  deformed  CFD  meshes  computed  by  the  BEM  solver  are  shown 
in  figure  6.7. 


6.2  Deforming  GFD  Mesh  of  the  ARW-2  Wing 

The  ARW-2  wing  (the  second  of  Aeroelastic  research  wing)  was  developed  by  NASA  Langley 
and  tested  in  the  transonic  dynamic  tunnel  (TDT)  in  the  1980’s.  The  Nastran  structural 
model  of  the  ARW-2  wing  is  converted  from  an  original  spar-based  FEM  model  reported  in 
Ref  [96]. 

Figure  6.2  presents  the  Nastran  FEM  of  the  ARW-2  wind  tunnel  model  whose  natural 
frequencies,  generalized  masses,  and  generalized  stiffness  are  shown  in  table  6.2. 

Shown  in  Figure  6.8  is  the  surface  CFD  mesh  of  the  ARW-2  wing  which  consists  of  153x53 
grip  points.  Similar  to  the  generation  of  defomed  sufrace  CFD  mesh  of  the  MAVRIC-I 
wing,  the  BEM  solver  is  adopted  to  compute  the  deformed  surface  CFD  mesh  of  the  ARW-2 
wing  for  the  first  eight  modes.  Figure  6.9  depicts  the  deformed  FEM  modes  computed  by 
MSC.Nastran  and  the  deformed  CFD  mesh  computed  by  the  BEM  solver.  It  can  be  seen 
that  these  two  sets  of  deformed  models  are  very  similar.  Also,  the  deformed  CFD  meshes 
appear  to  smooth;  showing  the  robustness  of  the  BEM  solver. 
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Figure  6.8:  Deforming  FEM  Models  and  CFD  Meshes  of  the  8  MAVRIC-I  Wing  Structural 


Chapter  7 

Duke  3D  Wings  LCO:  AGARD  445.6 


Summary 

•  The  inviscid  Euler  CFD/HB  solution  method  has  been  extended  to  three-dimensional 
flow  over  a  wing.  The  AGARD  445-6  has  been  used  as  an  example  to  demonstrate  the 
capability  of  the  method.  To  our  knowledge  these  are  the  first  LCO  results  obtained 
for  this  wing  by  any  investigators.  Among  the  principal  findings  is  that  overa  wide 
Mach  number  range  there  is  no  stable  LCO  that  would  be  detectable.  Moreover  in 
the  low  supersonic  range  an  unstable  LCO  is  predicted  which  suggests  that  LCO  of  a 
large  amplitude  would  be  encountered  below  the  predicted  flutter  boundary.  This  is  an 
interesting  result  because  previous  correlations  of  experiment  and  theory  have  suggested 
that  flutter  was  observed  experimentally  at  dynamic  pressures  below  the  predicted  flutter 
boundary  at  low  supersonic  Mach  numbers.  The  unstable  LCO  predicted  at  these  Mach 
number  may  be  at  least  part  of  the  explanation  of  this  difference  between  experiment 
and  linear  flutter  theory  which  has  been  so  perplexing  to  the  aerospace  community. 

•  The  complete  work  in  detail  can  be  found  in  the  paper  entitled,  ”A  Harmonic  Bal¬ 
ance  Approach  to  Modeling  Three-Dimensional  Nonlinear  Unsteady  Aerodynamics  and 
Aeroelasticity,”  by  J.P.  Thomas,  E.H.  Dowell,  and  K.C.  Hall  presented  at  the  ASME 
International  Mechanical  Engineering  Conference  and  Exposition,  ASME  Paper  IMECE- 
2002-32532,  November  2002  in  New  Orleans,  LA.  It’s  abstract  reads  as  follows: 
Presented  is  a  frequency  domain  harmonic  balance  (HB)  technique  for  modeling  non¬ 
linear  unsteady  aerodynamics  of  three-dimensional  transonic  inviscid  flows  about  wing 
configurations.  The  method  can  be  used  to  model  efficiently  nonlinear  unsteady  aero¬ 
dynamic  forces  due  to  finite  amplitude  motions  of  a  prescribed  unsteady  oscillation  fre¬ 
quency.  When  combined  with  a  suitable  structural  model,  aeroelastic  (fluid- structure), 
analyses  may  be  performed  at  a  greatly  reduced  cost  relative  to  time  marching  methods 
to  determine  the  limit  cycle  oscillations  (LCO)  that  may  arise.  As  a  demonstration 
of  the  method,  nonlinear  unsteady  aerodynamic  response  and  limit  cycle  oscillation 
trends  are  presented  for  the  AGARD  445-6  wing  configuration.  Computational  results 
based  on  the  inviscid  flow  model  indicate  that  the  AGARD  445-6  wing  configuration  ex¬ 
hibits  only  mildly  nonlinear  unsteady  aerodynamic  effects  for  relatively  large  amplitude 
motions.  Furthermore,  and  most  likely  a  consequence  of  the  observed  mild  nonlinear 
aerodynamic  behavior,  the  aeroelastic  limit  cycle  oscillation  amplitude  is  predicted  to 
increase  rapidly  for  reduced  velocities  beyond  the  flutter  boundary.  This  is  consistent 
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with  results  from  other  time-domain  calculations.  Although  not  a  configuration  that 
exhibits  strong  LCO  characteristics,  the  AGARD  445.6  wing  nonetheless  serves  as  an 
excellent  example  for  demonstrating  the  HB/LCO  solution  procedure. 

•  Presented  below  are  two  essential  conclusive  figures  showing  the  HB  flutter/LCO  solu¬ 
tions  due  to  our  inviscid  (Euler)  model  of  a  3D  AGARD  445-6  wing.  The  left  figure 
shows  that  the  HB  result  in  flutter  boundary  throughout  the  transonic  Mach  numbers  in 
which  it  correlates  well  with  test  data  in  the  lower  transonic  Mach  number  range.  The 
right  figure  shows  that  the  HB  method  predicts  and  unstable  LCO  in  the  low  supersonic 
Mach  number  range,  suggesting  the  LCO  of  a  large  amplitude  would  be  encountered 
below  the  predicted  flutter  boundary. 

7.1  Introduction 

Recently,  Hall  et.  al.  [78]  presented  a  new  and  and  computationally  efficient  harmonic 
balance  (HB)  approach  for  modeling  nonlinear  periodic  unsteady  aerodynamics  of  two- 
dimensional  turbo-machinery  configurations.  The  HB  method  can  be  used  to  model  non¬ 
linear  unsteady  aerodynamic  response  due  to  finite  amplitude  motions,  and  it  is  easily  im¬ 
plemented  within  the  framework  of  a  conventional  iterative  CFD  method.  Subsequently, 
the  present  authors  [?,  80]  devised  a  computational  strategy  based  on  the  harmonic  balance 
technique  for  modeling  aeroelastic  limit  cycle  oscillation  (LCO)  behavior  of  two-dimensional 
airfoil  configurations. 

In  the  following  article,  we  now  extend  both  the  harmonic  balance  technique  for  modeling 
nonlinear  unsteady  aerodynamics,  and  the  harmonic  balance  based  LCO  solution  procedure, 
to  three-dimensional  inviscid  flows  about  wing  configurations.  The  main  challenge  in  going 
from  two  to  three-dimensions  is  that  structural  models  in  three-dimensions  normally  consist 
of  many  more  degrees-of-freedom  (DOF)  other  than  the  two  “typical”  pitch  and  plunge  DOF 
for  two-dimensional  airfoil  sections. 

For  two-dimensional  flows,  it  has  been  observed  that  a  Newton-Raphson  root  finding  tech¬ 
nique  in  conjunction  with  the  harmonic  balance  method  proves  to  be  a  very  efficient  pro¬ 
cedure  for  computing  LCO  solutions  [?,  80].  However  for  the  Newton-Raphson  LCO  so¬ 
lution  procedure,  one  must  compute  gradients  (or  sensitivities)  of  the  nonlinear  unsteady 
aerodynamic  loading  with  respect  to  each  of  the  structural  degrees-of-freedom.  For  a  two- 
dimensional  airfoil  section,  this  requires  only  two  gradient  calculations. 

For  three-dimensional  configurations  however,  this  step  may  have  to  be  carried  out  for  many 
more  structural  degrees-of-freedom.  One  focus  of  the  research  presented  in  this  paper  has 
been  to  determine  to  what  extent  one  may  reduce  the  computational  effort  by  using  a  lim¬ 
ited  number  of  the  possible  structural  DOF.  As  will  be  shown  subsequently,  using  only  the 
structural  modes  that  dominate  the  flutter  onset  condition  appears  to  be  sufficient  to  model 
LCO  behavior  accurately  when  using  the  HB/LCO  solution  procedure. 
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Figure  7.1:  Mach  Number  Flutter  Speed  Dip  Trend  for  the  AGARD  445.6  Wing  Configura¬ 
tion 


In  the  following,  we  demonstrate  the  harmonic  balance  technique  for  modeling  nonlinear 
unsteady  aerodynamics,  along  with  the  harmonic  balance  based  LCO  solution  procedure, 
for  the  well  known  AGARD  445.6  transonic  wing  configuration  [?,  82].  We  begin  by  present¬ 
ing  the  aeroelastic  system  of  equations  applicable  to  the  AGARD  445.6  wing  configuration. 


7.2  Aeroelastic  Model  Governing  Equations 

The  governing  aeroelastic  system  of  equations  for  the  AGARD  445.6  wing  configuration 
(Fig.  7.3)  may  be  written  in  the  frequency  domain  as: 


-  Cq(£,w)  =0. 


For  the  present  analysis,  we  consider  a  linear  structural  model  expressed  in  terms  of  struc¬ 
tural  natural  modes  together  with  a  nonlinear  aerodynamic  model.  Structural  damping  is 
neglected.  kw  is  a  constant  dependent  on  the  wing  shape  and  overall  mass  given  by 


km  — 


7T  Ar  ( 1-t-A  t )  ( l-hAH~-A  t2 ) 


and  Cq  is  the  vector  of  generalized  aerodynamic  forces,  the  rrfi1  element  of  which  is  given  by 
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Figure  7.2:  Computed  AGARD  445.6  Wing  Configuration  LCO  Behavior  Trends 


Ce"'=^JLp^-&dA  (7-3) 

where  the  integral  is  evaluated  over  the  surface  of  the  wing. 

In  addition  to  the  generalized  mass  matrix  Al,  wing  structural  frequency  ratio  matrix  f2, 
and  mass  ratio  /x,  the  governing  aeroelastic  equations  also  include  the  reduced  velocity  V, 
reduced  frequency  u,  and  structural  modal  coordinates  £. 

As  was  done  for  the  aeroelastic  LCO  model  of  Ref.  [80],  one  can  also  include  the  static 
aeroelastic  equations  in  the  mathematical  model.  However  for  the  case  of  the  AGARD  445.6 
wing  configuration,  this  is  not  necessary  since  the  configuration  consists  of  a  wing  based  on 
a  constant  symmetric  airfoil  section  at  zero  degree  angle-of-attack.  As  such,  no  net  steady, 
or  static,  aerodynamic  load  acts  on  the  wing. 

As  will  be  discussed  in  the  following  section,  in  order  to  solve  the  aeroelastic  system  (Eq.  8.1), 
one  must  be  able  to  compute  the  generalized  aerodynamic  forces  C  q  for  finite  values  of  the 
structural  modal  coordinates  £.  This  is  where  the  harmonic  balance  procedure  becomes  an 
invaluable  tool. 
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7.3  Harmonic  Balance  Method 

7.3.1  Governing  Equations 

We  consider  the  inviscid  Euler  equations  (the  Reynolds  averaged  Navier-Stokes  equations 
can  be  treated  in  a  similar  manner),  which  may  be  written  in  integral  form  as 


d_ 

dt. 


/fcv+/i^-Ux')-fidA=° 


where  U  is  the  vector  of  conservative  fluid  variables 


U  =  {p  pu  pv  pw  Et}2 


and 


T  =  Ft  +  Gj  +  nk 

where  F,  G,  and  H  are  the  x,  y,  and  2:  direction  component  flux  vectors,  i.e. 

F= 


/  \ 
pu 

f 

pv 

f  \ 

pw 

pu2jt-p 

puv 

puw 

puv 

►  ,  G=  < 

pv2+p 

>,  H=< 

puv  ^ 

puw 

pvw 

pv2+p 

!Et+p)u 

iEt+p)v4 

K(Et+p)wj 

The  unsteady  motion  of  the  control  volume  x  is  given  by 

x  =  fi  +  gj  +  hk, 


(7.4) 


(7.5) 


(7.6) 


(7.7) 


and  this  accounts  for  the  effect  of  wing  surface  and  grid  motion. 


7.3.2  Fourier  Series  Expansion 

We  consider  the  unsteadiness  of  the  flow  to  be  strictly  periodic  in  time  with  period  T  =  2tt/u> 
where  u  is  the  fundamental  unsteady  frequency.  Thus  we  can  expand  Eq.  7.4  in  a  Fourier 
series.  For  example, 


III; 


Nh 


UdV  =  Q(t)«  V  Q„^ 
v(t)  u=_Nh 


(7.8) 


so  that 
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(7.9) 


d_ 

at 


Ilkr 


nh 

U=—Nh 


and  similarly 


f  f  (.F-Ux)  -iidA  =  R(t)  «  V  ’RnejrUjJt. 


Nh  is  the  number  of  harmonics  used  in  the  Fourier  expansion. 


(7.10) 


7.3.3  Fourier  Coefficients 

Substituting  the  Fourier  expansions  (Eqs.  7.9  and  7.10)  into  Eq.  7.4,  multiplying  by  e-jrrujt, 
and  integrating  over  one  period,  for  each  m,  i.e. 

/  m  E  (jnQn  +  BjjeP^e-^dt  (7.11) 

J°  1  n=—Njj  J 

yields  a  system  of  equations  for  the  Fourier  coefficients.  Namely 

AQ  +  R  =  0  (7.12) 

where 


A  = 

-jNH 

,Q=  < 

Q -nh  ' 

Q-Nh+1 

>  ,R=  < 

'  R -nh  ' 
R-Wtf+i 

3nh 

.  Q  Nh  j 

.  R  nh  j 

7.3.4  Time  Domain  Variables 

Next,  via  a  Fourier  transform  matrix  E,  one  can  relate  the  Fourier  coefficient  variables  Q 
to  time  domain  solution  variables  stored  at  uniformly  space  sub-time  levels  within  a  given 
period  of  motion,  i.e., 


Q  —  EQ,  R  —  ER 


(7.13) 


where 
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) 


(7.14) 


and 


More  specifically, 


and 


Thus 


Q  =  < 

Q(t0)  } 
Q(ti) 

> ,  R  =  < 

r  R(t0)  ' 
R(ti) 

kQ(^2ATH)  j 

R(t2NH)j 

27rn 

(2  Njj  +  l)u 


n 


0, 1,  ...,2Nh- 


'  fffv(t0)U(to)dV  ' 

///„<«, >U(«dV 

Q=  <  .  >, 

ffA(t0)  f^(to)-U(to)x(t0))  n(t0)dA  ' 

(■^■(*i)-U(ti)x(ti)J-n(ti)dA 

.  If  A(t2 Nh  )  {£(t2NH  )  -  u (t2NH  )x(t2NH  ))  •  n(i2ATH )  dA 


(7.15) 


(7.16) 


(7.17)  ' 


AEQ  -)-  ER  —  0, 


(7.18) 


and 


E-1AEQ  +  E_1ER  =  0. 


(7.19) 


So  now  one  can  work  in  terms  of  the  time  domain  variables,  which  is  in  general  much  easier 
to  do.  Note  however  that  we  take  full  advantage  of  the  assumed  periodic  motion  and  thus 
a  transient  time  simulation  and  its  associated  computational  cost  is  avoided.  The  resulting 
system  of  equations  can  then  be  written  as 

DQ  +  R  =  0  (7.20) 

where 

D  =  E_1AE.  (7.21) 
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7.3.5  Pseudo  Time  Marching 

By  adding  a  pseudo  time  derivative  term  5Q/St  to  Eq.  7.20,  one  can  then  develop  an  itera¬ 
tive  technique  for  determining  the  harmonic  balance  solution  Q.  That  is, 

^  +  DQ  +  R  =  0,  (7.22) 

whereby  one  then  simply  “marches”  Eq.  7.22  in  a  fictitious  time-like  manner  until  a  steady 
state  (Eq.  7.20)  is  achieved.  Solution  acceleration  techniques  such  as  local  time-stepping, 
pre-conditioning,  residual  smoothing,  and  multi-grid  can  all  be  used  to  accelerate  the  con¬ 
vergence  of  the  harmonic  balance  solution. 

So  for  example,  in  the  case  of  a  finite-volume  based  CFD  method,  Eq.  7.22  is  solved  for 
every  computational  “finite-volume”  comprising  the  computational  mesh.  The  overall  thus 
consists  of  pseudo  time  marching  N  x  (27V#  + 1)  dependent  variables  where  N  is  the  number 
of  mesh  points  times  the  number  of  dependent  variables. 

Modifying  an  existing  steady  CFD  flow  solver  to  implement  the  harmonic  balance  technique 
is  thus  a  relatively  straight  forward  task  as  the  main  requirement  is  just  a  re-dimensioning 
of  the  primary  arrays  from  N  elements  to  N  x  (2 NH  +  1)  elements.  Again  since  solution 
acceleration  techniques  such  as  local  time-stepping,  pre-conditioning,  residual  smoothing, 
and  multi-grid  can  be  used,  the  computational  cost  of  the  method  in  determining  unsteady 
solutions  thus  scales  by  a  factor  of  (2Nh+1)  times  the  cost  of  a  nominal  steady  flow  solution. 

Once  a  harmonic  balance  solution  has  been  determined  for  a  pre-specified  wing  motion  and 
frequency,  one  can  then  obtain  the  resulting  generalized  aerodynamic  forces  from  Eq.  8.3. 


7.4  Limit  Cycle  Oscillation  Solution  Procedure 

Based  on  previous  research  efforts  conducted  for  two-dimensional  airfoil  configurations,  we 
have  found  that  a  Newton-Raphson  technique  in  conjunction  with  the  harmonic  balance 
nonlinear  unsteady  flow  solver  technique  provides  an  efficient  method  to  determine  limit 
cycle  oscillation  response. 

As  demonstrated  in  Refs.  [79]  and  [80],  the  harmonic  balance  LCO  solution  method  pro¬ 
ceeds  by  choosing  one  of  the  structural  modal  coordinates  to  be  the  independent  variable. 
For  the  two-dimensional  airfoil  case,  this  was  chosen  to  be  the  pitch  coordinate  a.  For 
three-dimensional  flow  about  wings,  we  now  consider  £i,  the  modal  coordinate  of  the  first 
(typically  bending)  structural  mode  shape,  as  the  independent  variable.  We  also  chose  £i  to 
be  real  valued. 

In  formulating  the  HB/LCO  solution  technique,  one  then  proceeds  by  dividing  Eq.  8.1 
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through  by  £1,  and  re-expressing  the  system  of  equations  as: 


R(L,6)=fc.A<(-oVl+  p«)  |  -  Cc(|-,6,0)  =  0.  (7.23) 

Considering  both  the  real  and  imaginary  parts  of  Eq.  7.23,  the  vector  L  then  represents  the 
unknown  LCO  solution  variables 

V 

u 

Re(6)/£i 

M6)/&  ,  (7.24) 

which  consists  of  the  LCO  reduced  velocity  V  (this  variable  is  customarily  the  independent 
variable  in  time-domain  LCO  solution  techniques),  LCO  reduced  frequency  u,  and  the  real 
and  imaginary  parts  of  the  ratio  of  each  of  the  LCO  structural  modal  coordinates,  for  modes 
two  and  higher,  to  the  first  structural  modal  coordinate  amplitude  £i.  The  real  and  imag¬ 
inary  parts  of  Eq.  7.23  as  such  represent  a  system  of  2 M  equations  for  the  2 M  unknown 
LCO  solution  variables  of  L. 


As  observed  in  Refs.  [79]  and  [80],  an  efficient  method  for  solving  Eq.  7.23  is  to  use  a 
simple  Newton-Raphson  root  finding  technique.  This  is  an  iterative  method  for  solving  for 
the  unknown  LCO  variables  L  whereby  one  “marches”  the  vector  equation 

Ln+1  =  Ln  —  R(Ln),  (7.25) 

until  a  suitable  level  of  convergence  is  achieved. 


We  have  also  observed  that  one  can  use  simple  forward  finite-differencing  to  compute  the 
column  vectors  of  <9R(L) /dL.  That  is, 


<9R(L)~ 

dL 


SR  SR  SR _ dR 

dv  d“  aR e(£j)/£i  almfe)/?i 


(7.26) 


where  for  example 


9R(L) 

dV 


R(L,  V  +  e)  -  R(L,  V) 
e  ' 


100 


etc.  for  a  small  e.  Determining  the  column  vectors  of  <9R/<9L  in  this  manner  thus  requires 
numerous  computations  of  R(L)  for  various  perturbations  to  L.  This  in  turn  means  several 
computations  of  the  unsteady  aerodynamic  loading  C  q  for  the  different  perturbations  of 
L,  and  this  is  where  the  harmonic  balance  solver  comes  into  play.  This  is  also  the  most 
computationally  expensive  aspect  of  the  LCO  solution  methodology.  However,  the  overall 
method  does  lend  itself  to  a  simple  computational  parallelization  strategy  as  each  of  the 
gradient  approximations  can  be  calculated  on  an  individual  computer  processor. 

7.4.1  The  AGARD  445.6  Transonic  Wing  Configuration 

As  noted  in  the  introduction,  in  order  to  demonstrate  the  harmonic  balance  nonlinear  aero¬ 
dynamic  and  LCO  solution  procedures,  the  the  AGARD  model  445.6  transonic  wing  con¬ 
figuration  [?,  82]  is  chosen  as  an  example.  This  is  a  45  degree  quarter  chord  swept  wing 
based  on  the  NACA  64A004  airfoil  section,  which  has  an  aspect  ratio  of  3.3  (for  the  full 
span),  and  a  taper  ratio  of  2/3.  Figure  7.3  illustrates  the  computational  mesh  employed 
for  this  configuration  with  Fig.  7.3a  showing  a  close-up  of  the  wing  surface  and  symmetry 
boundary  grids,  and  Fig.  7.3b  showing  the  outer  boundary  grid.  The  grid,  which  in  this 
instance  is  the  “medium”  resolution  mesh,  is  based  on  an  “O-O”  topology  that  employs  49 
(mesh  i  coordinate)  computational  nodes  about  the  wing  in  the  stream-wise  direction,  33 
(mesh  j  coordinate)  nodes  normal  to  the  wing,  and  33  (mesh  k  coordinate)  nodes  along  the 
semi-span.  The  total  number  of  fluid  dynamic  DOF  for  this  CFD  mesh  is  thus  266,805  (i.e. 
5dependent  flow  variables  x  imax  (49)  x  jmax( 33)  x  fcmax( 33)).  The  outer  boundary  of  the 
grid  extends  five  semi-spans  from  the  mid-chord  of  the  wing  at  the  symmetry  plane.  The 
particular  structural  configuration  of  the  wing  under  consideration  is  referred  to  as  the  “2.5 
ft.  weakened  model  3”  [?,  82]. 

Finally,  Fig  7.4  shows  the  first  three  computed  (via  a  finite  element  analysis)  structural 
mode  shapes  and  natural  frequencies  for  the  AGARD  445.6  “weakened”  wing  configuration 
as  presented  in  [81].  The  first  mode  shape  is  seen  to  consist  of  a  first  bending  type  of  mo¬ 
tion,  the  second  mode  shape  then  being  a  first  twisting  type  of  motion,  and  finally  the  third 
mode  being  a  second  bending  motion.  Two  additional  mode  shapes,  second  twisting  and 
third  bending,  are  also  provided  in  [81].  These  mode  shapes  have  all  been  mapped  to  the 
computational  meshes  used  in  this  study  using  a  sixth-order  multi-dimensional  least  squares 
curve  fitting  technique. 


7.4.2  Unsteady  Grid  Motion  Treatment 

As  noted  in  section  3.4,  for  the  harmonic  balance  method,  the  flow  solution  variables  are 
stored  at  2NH  +  1  sub-time  levels  over  a  period  of  one  cycle  of  motion.  This  means  that 
the  unsteady  deformed  shape  of  the  wing  is  thus  required  at  2 Nh  +  1  sub-time  levels  over 


(a)  Wing  Surface  and  Symmetry  Plane  Grids 


(b)  Far- field  Boundary  Grids 


Figure  7.3:  Computational  Grid  Layout  for  AGARD  445.6  Wing  Configuration.  49(i)  x  33(j) 
x  33  (k)  Grid  Shown. 
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Structural  Mode  One  ( First  Bending)  -  60  hz 


Structural  Mode  Three  (Second  Bending)  -  304  Hz 
Figure  7.4:  AGARD  445.6  Wing  “Weakened”  Structural  Configuration  Mode  Shapes. 


the  same  period  of  motion.  In  the  following  analysis,  we  consider  a  linear  structural  model 
whereby  we  approximate  the  finite-amplitude  motions  as  a  superposition  of  a  limited  number 
of  the  natural  mode  shapes  of  the  structure.  That  is, 

(7.27) 

where  x0.  j  fc  is  the  nominal  stationary  grid  and  (f>i}k(j)  is  a  blending  function,  which  is  equal 
to  one  at  the  wing  surface  boundary,  and  which  decays  to  zero  at  the  outer  far-field  grid 
boundary.  In  this  instance,  we  use  the  following  blending  function 

=  1  —  Si,k(j)  /  Si,k{jmax)  (7.28) 

where  si)k(j)  is  the  distance,  in  the  j^1  mesh  coordinate,  along  a  grid  line  curve  emanating 
from  the  (i,  fc)^1  mesh  point  on  the  surface  of  the  wing. 

This  procedure  for  modeling  the  finite  amplitude  wing  surface  and  grid  motion  is  very  simple 
and  efficient.  However,  it  can  eventually  lead  to  negative  volumes  when  the  amplitude  of  the 
motion  is  too  great.  This  is  in  part  due  to  the  “O-O”  grid  topology,  which  has  the  tendency 
to  have  a  rather  highly  skewed  mesh  near  the  wing  tip. 


M 
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7.5  Nonlinear  Unsteady  Aerodynamics 

We  next  compute  the  nonlinear  aerodynamic  response  for  a  Moo=0.960  background  steady 
flow  about  the  AGARD  445.6  wing  configuration  due  to  various  finite  amplitudes  of  unsteady 
wing  motion.  We  choose  a  reduced  frequency  of  u  =  0.1  since  this  is  near  the  reduced  fre¬ 
quency  at  which  the  wing  flutters. 


7.5.1  Mesh  Convergence 

First,  we  consider  the  quality  of  the  mesh.  To  get  a  sense  of  required  mesh  resolution, 
we  consider  a  small,  in  effect  linear,  amplitude  motion  of  the  first  modal  coordinate,  i.e. 

=  0.0001  for  a  reduced  frequency  CD  =  0.1.  Of  course  grid  resolution  is  also  an  impor¬ 
tant  issue  for  larger  amplitude  motions.  However  for  larger  motion  amplitudes,  the  number 
of  harmonics  one  uses  in  the  HB  method,  in  addition  to  the  grid  resolution,  plays  a  role  in 
overall  solution  accuracy.  As  such,  a  comprehensive  study  of  model  resolution  becomes  quite 
a  bit  more  involved.  For  this  reason,  we  have  chosen  to  examine  a  small  amplitude  motion. 
Future  grid  convergence  studies  will  consider  larger  amplitude  motions  and  higher  harmonics. 

Figure  7.5  shows  the  computed  steady  and  unsteady  flow  surface  pressure  distributions  at 
70%  span  for  three  different  grid  resolutions.  As  can  be  seen,  there  is  not  much  difference  in 
the  solutions  amongst  the  three  different  grid  resolutions.  As  Such,  we  have  chosen  to  employ 
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Mean  ( Zeroth  Harmonic)  Surface  Pressure 


Real  Unsteady  (First  Harmonic)  Surface  Pressure 


Imaginary  Unsteady  (First  Harmonic)  Surface  Pressure 


Figure  7.5:  Mesh  Convergence  Characteristics  for  the  AGARD  445.6  Wing  Surface  Pressure 
at  70%  of  Wing  Semi-span:  a0  =  0.0  (deg),  Moo  =  0.960,  Hi  =  0.1,  and  &  =  0.0001. 
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Structural  Mode  Shape 

zWv/(&/2) 

^imax  /  °r 

1  -  First  Bending 

0.2416 

IHsEEll 

2  -  First  Twisting 

0.0500 

2  -  Second  Bending 

0.0030 

0.0547 

0.0749 

Table  7.1:  AGARD  445.6  Wing  Configuration  Maximum  Modal  Coordinate  Amplitudes  - 
49  (i)  x  33  (j)  x  33  (k)  Grid. 


49x33x33  grid  for  all  subsequent  calculations.  As  can  be  seen,  since  the  AGARD  445.6  wing 
uses  only  a  4%  thick  symmetric  airfoil  section,  the  pressure  distribution  is  also  symmetric 
about  the  upper  and  lower  surfaces,  and  because  the  section  is  so  thin,  no  definitive  shock 
is  readily  apparent  at  even  this  high  transonic  Mach  number. 

7.5.2  Finite  Amplitude  Motion  Limits 

As  mentioned  in  section  4.2,  negative  computational  cell  volumes  can  occur  once  the  ampli¬ 
tude  of  the  wing  motion  becomes  too  large.  Table  7.1  shows  the  maximum  value  for  each  of 
the  first  three  structural  modal  coordinates,  when  taken  separately,  at  which  point  negative 
volume  problems  start  to  occur  for  the  49x33x33  mesh.  Also  listed  is  the  corresponding 
maximum  unsteady  vertical  displacement  of  the  wing  surface  zimax  during  an  interval  of 
motion. 


7.5.3  Harmonic  Convergence 

Next,  Figs.  7.6  and  7.7  show  computed  real  and  imaginary  parts,  respectively,  of  the  first 
three  generalized  force  coefficients  (Cqv  Cq2,  and  ca3)  normalized  by  the  each  modal  coor¬ 
dinate  (£i,  £2,  and  £3)  as  a  function  of  the  magnitude  of  each  modal  coordinate  (£1,  £2,  and 
£3).  Shown  are  results  when  using  one,  two,  and  three  harmonics  in  the  harmonic  balance 
solution  procedure.  As  can  be  seen,  the  amplitude  of  the  modal  coordinate  has  an  effect  on 
the  generalized  force,  and  particularly  so  for  the  real  part.  Also  evident  is  how  one  should  use 
more  harmonics  when  considering  ever  increasing  modal  coordinate  amplitudes.  However, 
in  this  example,  the  use  of  only  two  harmonics  can  be  seen  to  produce  a  very  high  level  of 
harmonic  convergence  for  even  the  largest  modal  amplitudes.  It  has  been  our  experience 
that  two  harmonics  in  many  cases  are  all  that  are  necessary  for  achieving  a  good  level  of 
accuracy  as  long  as  the  amplitudes  are  not  too  large,  say  for  example,  a  maximum  of  five 
degrees  in  pitch  amplitude  for  the  case  of  an  airfoil  section. 

Next,  Fig.  7.8  shows  surface  pressure  distributions  for  the  zeroth  harmonic  (Fig.  7.8a)  (i.e. 
mean  flow),  along  with  the  real  (Fig.  7.8b)  and  imaginary  (Fig.  7.8c)  parts  of  the  first  har¬ 
monic  (unsteady  flow),  again  normalized  by  the  modal  coordinate  amplitude  £1  at  a  location 
approximately  70%  of  wing  semi-span.  In  this  instance,  we  consider  four  different  modal 
coordinate  amplitudes  (£i=0.0001,  0.0020,  0.0050,  0.0100)  for  the  case  of  three  harmonics 
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Structural  Mode  One 


Structural  Mode  Two 


Structural  Mode  Three 


Figure  7.6:  Real  part  of  Unsteady  Generalized  Force  as  a  Function  of  Mode  Shape  Amplitude: 
<*o  =  0.0  (deg),  Moo  =  0.960,  and  u>  =  0.1. 
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Structural  Mode  One 


Structural  Mode  Two 


Structural  Mode  Three 


Figure  7.7:  Imaginary  part  of  Unsteady  Generalized  Force  as  a  Function  of  Mode  Shape 
Amplitude:  a0  =  0.0  (deg),  =  0.960,  and  u>  =  0.1. 
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Mean  (Zeroth  Harmonic)  Surface  Pressure 


Real  Unsteady  (First  Harmonic)  Surface  Pressure 


Imaginary  Unsteady  (First  Harmonic)  Surface  Pressure 

Figure  7.8:  Mean  and  Unsteady  Flow  Surface  Pressure  Distribution  Dependence  on  Mode 
Shape  Amplitude:  a0  =  0.0  (deg),  Moo  =  0.960,  Q  =  0.1,  and  NH  =  3. 
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used  in  the  harmonic  balance  technique.  Readily  apparent  are  nonlinear  effects  due  to  the 
ever  increasing  modal  coordinate  amplitudes.  In  fact, this  is  surprisingly  so  for  the  zeroth 
harmonic  of  the  solution  (Fig.  7.8a). 

Finally,  Fig.  7.9  shows  nor realized  zeroth  and  first  harmonic  surface  pressure  distributions 
for  the  largest  modal  coordinate  amplitude.  (fi=0.0120)  while  varying  the  number  of  har¬ 
monics  used  in  the  harmonic  .balance  technique.  Two  harmonics  can  be  seen  to  produce 
good  harmonic  convergence. 

7.6  Nonlinear  Unsteady  Aeroelasticity 

7-6-1  Flutter  Onset  Trend  for  the  AG ARD  445.6  Wing  Configu¬ 
ration 

Figure  7.-10;  shows  the  •i^ffld>^ed*Mach-'>n'Ufidljer  flutter  onset  reduced  velocity  trend  fOi  the 
AGARD  445.6  wing  configuration.  Computational  model  results  are  shown  together;  with 
experimental  results.  The  six  Mach  numbers  are  the  same  as  those  of  the  experimental  in¬ 
vestigation  of  Ref.  [82].  Tfee^reduction-  of  the  flutter  onset  velocity  in  the  transonic  region  is 
readily  apparent.  The  flutter  onset  conditions  are  determined  using  the  proper  orthogonal 
decomposition  (POD)  reduced  order  model  (ROM)  time-linearized  unsteady  aerodynamic 
approach  of  Ref.  [84].  The*  POD /ROM  approach  has  been  demonstrated  to  provide  a  good 
approximation  of  the  flutter  onset  conditions,  which  can  then  be  used  as  initial  conditions 
for  the  iterative  HB/LCO  solution  procedure.  The  agreement  between  theory  and  experi¬ 
ment  is  generally  good  except  at  the  supersonic  free-stream  Mach  numbers.  Other  imdscid 
computational  studies  by  Investigators  such  as  Lee-Rausch  and  Batina  [85]  and  Gordnier 
and  Melville  [86]  have  demonstrated  similar  differences  between  theory  and  experiment  for 
the  supersonic  Mach  numbers. 

Table  7.2  shows  computed  Vhlues  for  the  reduced  velocity  V,  reduced  frequency  d>,  and  the 
shape  of  the  unsteady  motion  of  the  wing-at  the  flutter  onset  condition,  which  is  presented  in 
terms  of  the  ratio  of  the  amplitude  of  each  of  structural  modal  coordinates  normalized  to  the 
amplitude  of  first  structural 'modal  coordinate.  As  can  be  seen,  the  first  bending  structural 
motion  dominates  the  flutter  -onset  unsteady  motion  for  the  six  Mach  numbers  considered, 
particularly  in  the  transoniefregdOn.  In  faict;  beyond  the  first  three  structural  modes,  higher 
modes  Only  contribute  a  fraction' of  a  percent  to  the  overall  unsteady  motion  at  the  flutter 
onset  condition.  As  will  be 'Shown,  one  can  typically  neglect  these  higher  numbered  modes 
in  the  HB/LCO  solution  procedure. 


7.6.2  LCO  Behavior  itends  for  the  AGARD  445.6  Wing  Config¬ 
uration 


no 


Airfoil  Surface  Location,  x/c 


Mean  (Zeroth  Harmonic)  Surface  Pressure 


Real  Unsteady  (First  Harmonic)  Surface  Pressure 


Airfoil  Surface  Location,  x/c 


Imaginary  Unsteady  (First  Harmonic)  Surface  Pressure 

Figure  7.8:  Mean  and  Unsteady  Flow  Surface  Pressure  Distribution  Dependence  on  Mode 
Shape  Amplitude:  a0  =  0.0  (deg),  =  0.960,  a>  =  0.1,  and  NH  =  3.  ' 
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used  in  the  harmonic  balance  technique.  Readily  apparent  are  nonlinear  effects  due  to  the 
ever  increasing  modal  coordinate  amplitudes.  In  fact, this  is  surprisingly  so  for  the  zeroth 
harmonic  of  the  solution  (Fig.  7.8a).  ; 

Finally,  Fig.  7.9  shows  normalized  zeroth  and  first  harmonic  surface  pressure  distributions 
for  the  largest  modal  coordinate  amplitude  (£i=0.0120)  while  varying  the  number  of  har¬ 
monics  used  in  the  harmonic  balance  technique.  Two  harmonics  can  be  seen  to  produce 
good  harmonic  convergence. 


7.6  Nonlinear  Unsteady  Aeroelasticity 

7.6,1  Flutter  Onset  Trend  for  the  AGARD  445,6  Wing  Configu¬ 
ration 

Figure  7.10  shows  the  computed  Mach  number  flutter  onset  reduced  velocity  trend  for  the 
AGARD  445.6  wing  configuration.  Computational  model  results  are  shown  together  with 
experimental  results.  The  six  Mach  numbers  are  the  same  as  those  of  the  experimental  in¬ 
vestigation  of  Ref.  [82].  The  reduction  of  the  flutter  onset  velocity  in  the  transonic  region  is 
readily  apparent.  The  flutter  onset  conditions  are  determined  using  the  proper  orthogonal 
decomposition  (POD)  reduced  order  model  (ROM)  time-linearized  unsteady  aerodynamic 
approach  of  Ref.  [84].  The  POD/ROM  approach  has  been  demonstrated  to  provide  a  good 
approximation  of  the  flutter  onset  conditions,  which  can  then  be  used  as  initial  conditions 
for  the  iterative  HB/LCO  solution  procedure.  The  agreement  between  theory  and  experi¬ 
ment  is  generally  good  except  at  the  supersonic  free-stream  Mach  numbers.  Other  inviscid 
computational  studies  by  investigators  such  as  Lee-Rausch  and  Batina  [85]  and  Gordnier 
and  Melville  [86]  have  demonstrated  similar  differences  between  theory  and  experiment  for 
the  supersonic  Mach  numbers. 

Table  7.2  shows  computed  values  for  the  reduced  velocity  V ,  reduced  frequency  u>,  and  the 
shape  of  the  unsteady  motion  of  the  wing  at  the  flutter  onset  condition,  which  is  presented  in 
terms  of  the  ratio  of  the  amplitude  of  each  of  structural  modal  coordinates  normalized  to  the 
amplitude  of  first  structural  modal  coordinate.  As  can  be  seen,  the  first  bending  structural 
motion  dominates  the  flutter  onset  unsteady  motion  for  the  six  Mach  numbers  considered, 
particularly  in  the  transonic  region.  In  fact,  beyond  the  first  three  structural  modes,  higher 
modes  only  contribute  a  fraction  of  a  percent  to  the  overall  unsteady  motion  at  the  flutter 
onset  condition.  As  will  be  shown,  one  can  typically  neglect  these  higher  numbered  modes 
in  the  HB/LCO  solution  procedure. 


7.6.2  LCu  Behavior  Trends  for  the  Agaru  445.6  Wing  Config¬ 
uration 


110 


Mach  Number,  M0 0 

0.499 

0.678 

0.901 

0.960 

1.072 

1.141 

v 

■KSSI 

0.450 

HKSSil 

miMii 

0.664 

u 

I 

0.141 

■KSBI 

0.0790 

nefcO/6 

-0.243 

-0.115 

-0.102 

-0.171 

0.0624 

0.0253 

■ 

0.0142 

Re(6)/6 

0.0492 

0.0360 

0.0179 

0.0905 

Im(6)/£i 

0.00379 

0.00202 

0.000389 

-0.000558 

-0.00246 

-0.00591 

Re(6i)/£i 

-0.00557 

-0.00477 

-0.00161 

-0.000110 

0.00707 

0.0117 

MU  )/U 

0.00495 

0.00346 

0.00203 

0.00104 

-0.00103 

-0.00166 

Re(&)/£i 

0.00102 

0.00101 

0.00117 

0.00116 

0.000437 

-0.00361 

Im(C5)/€i 

-0.00127 

-0.00067 

-0.000368 

-0.000486 

-0.000498 

-0.000574 

Table  7.2:  AGARD  445.6  Wing  Configuration  Mach  Number  Flutter  Onset  Conditions. 


Ill 


Mean  ( Zeroth  Harmonic)  Surface  Pressure 


Real  Unsteady  (First  Harmonic)  Surface  Pressure 


Imaginary  Unsteady  (First  Harmonic)  Surface  Pressure 


Figure  7.9:  Mean  and  Unsteady  Flow  Surface  Pressure  Distribution  Dependence  on  Number 
of  Harmonic  Used  in  HB  Solver:  a0  =  0.0  (deg),  =  0.960,  u>  =  o.l,  and  =  0.0120. 
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Figure  7.10:  Mach  Number  Flutter  Speed  Dip  Trend  for  the  AGARD  445.6  Wing  Configu¬ 
ration 


4JJ> 


Figure  7.11:  Computed  AGARD  445.6  Wing  Configuration  LCO  Behavior  Trends. 
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Figure  7.12:  Computed  AGARD  445.6  Wing  Configuration  LCO  Behavior  Trends  When 
Using  Different  Numbers  of  Structural  Modes  in  HB/LCO  Solutions  Procedure. 


Based  on  the  LCO  solution  procedure  presented  in  the  Section  4,  Fig.  7.11  shows  com¬ 
puted  LCO  behavior  trends  for  the  AGARD  445.6  wing  configuration  for  the  six  different 
Mach  numbers.  Unfortunately,  Ref.  [82]  does  not  provide  any  specific  experimental  LCO 
data.  Plotted  in  Fig.  7.11  is  the  LCO  amplitude  of  the  first  structural  modal  coordinate,  £1? 
versus  the  reduced  velocity  V  for  each  of  the  six  different  Mach  numbers.  The  data  points  for 
what  appear  to  be  zero  LCO  amplitude  actually  correspond  to  a  very  small  LCO  amplitude 
of  6  =  0.0001. 

In  this  instance,  only  the  first  three  structural  modes  are  used  for  the  HB/LCO  solution 
process.  The  convergence  of  the  HB/LCO  solution  procedure  is  also  very  rapid.  Only  one 
or  two  iterations  of  Eq.  7.25  are  required.  The  LCO  amplitudes  in  the  transonic  region  go 
to  larger  values  than  those  in  the  subsonic  and  high  supersonic  region  since  the  first  bending 
mode  shape  very  much  dominates  in  the  transonic  region.  As  such,  negative  volume  prob¬ 
lems  with  the  unsteady  grid  are  less  of  an  issue  in  the  transonic  region. 

Note  how  the  LCO  behavior  trend  curves  are  nearly  vertical.  This  is  indicative  of  very  linear 
aeroelastic  LCO  response  behavior,  which  is  not  all  that  surprising  for  this  particularly  thin 
wing  that  has  a  maximum  thickness  ratio  of  4%.  A  recent  time-domain,  and  viscous  flow, 
computational  study  of  this  same  configuration  for  Moo=0.960  and  Mo^l.141  conducted 
by  Gordnier  and  Melville  [86]  also  determined  that  no  LCO  conditions  appear  to  exist.  Note 
that  =  0.012  corresponds  to  an  unsteady  wing  amplitude  of  approximately  one  third  of 
the  wing  root  chord  (Section  5.2). 

Finally,  Fig.  7.12  shows  computed  LCO  behavior  trends  for  M ^=0.499  and  Moo=0.960 
when  using  two,  three,  and  four  structural  mode  shapes  in  the  HB/LCO  solution  procedure. 
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As  can  be  seen,  only  the  first  few  structural  mode  shapes  are  necessary  to  produce  converged 
LCO  solutions.  Again,  this  is  most  likely  due  to  the  fact  that  the  first  bending  mode  tends 
to  dominate  the  aeroelastic  flutter  motion  for  the  AGARD  445.6  wing. 

It  is  interesting  to  note  that  the  strongest  nonlinear  LCO  behavior  occurs  for  M00=1.Q72 
and  M00= 1.141  as  shown  in  Fig.  7.11.  For  these  Mach  numbers,  the  nonlinear  effect  is  to 
decrease  the  reduced  velocity  at  which  LCO  may  occur  below  the  flutter  onset  condition. 
This  may  be  at  least  a  partial  explanation  for  the  observed  differences  between  theory  and 
experiment  in  Fig.  7.10.  There  the  calculated  (linear)  reduced  velocity  at  the  onset  of  flut¬ 
ter  is  compared  to  the  reduced  velocity  at  which  flutter  was  observed  experimentally.  As 
shown  in  Fig.  7.10,  the  flutter/LCO  observed  experimentally  does  indeed  occur  at  a  reduced 
velocity  below  that  predicted  theoretically  for  the  onset  of  flutter. 


7.7  Conclusions 

A  harmonic  balance  method  for  modeling  nonlinear  periodic  unsteady  three-dimensional 
inviscid  transonic  flows  about  wing  configurations  is  presented.  Demonstrated  is  the  ability 
of  the  method  to  model  nonlinear  aerodynamic  effects  due  to  large  amplitude  motions  for 
a  Well  known  aeroelastic  configuration.  The  Use  of  only  two  harmonics  in  the  procedure  is 
shown  to  be  quite  sufficient  for  producing  harmonic  convergence.  A  limit  cycle  oscillation 
solution  methodology  is  also  presented  and  demonstrated  for  the  same  benchmark  transonic 
configuration.  A  weak  LCO  is  observed,  which  is  consistent  with  results  from  time-domain 
calculations  by  Gordnier  and  Melville  [86].  Nevertheless,  it  is  suggested  that  nonlinear 
aerodynamic  effects  may  at  least  partially  explain  some  of  the  previously  observed  differences 
between  experiment  and  linear  aeroelastic  theory. 


Nomenclature 


ar 

6,  c 

CQ 

Et 

f,9,h 

I 

j_ 

rh 

Moo 

M 

M 

Nh 


aspect  ratio  =  wing  span  squared/ wing  area 
semi-chord  and  chord  respectively 
vector  of  non-dimensional  generalized  forces 
total  energy 

scalar  functions  defining  unsteady  grid 
motion  in  x,  y,  and  2  coordinates,  respectively 
identity  matrix 

mass  of  wing 
free-stream  Mach  number 
number  of  structural  modes 
generalized  mass  matrix 

number  of  harmonics  used  in  harmonic  balance 
method 
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Pi 

Qoo 
P 5  Poo 

Q 

u,v,w 


Un 


V 


V 

<*0 

P 

$m_ 
M Ct 

n 


£nu  £ 


=  first  harmonic  unsteady  pressure 
=  free-stream  dynamic  pressure 
=  local  and  free-stream  density,  respectively 
=  vector  of  generalized  aerodynamic  forces 
=  flow  velocity  components  in  x,  y,  and  z  coordinate 
directions,  respectively 
=  free-stream  velocity 

=  volume  of  a  truncated  cone  having  stream-wise  root 
chord  as  lower  base  diameter,  stream- wise  tip  chord 
as  upper  base  diameter,  and  wing  half  span  as  height 
=  reduced  velocity  or  speed  index  V  =  Uoo/y/pVab 
=  airfoil  or  wing  root  steady  flow  angle-of-attack 
=  wing  taper  ratio  Xt  =  Ct/Cr 
=  mass  ratio,  n  =  m/p^v 
=  mfi1  structural  mode  shape 
=  frequency  and  reduced  frequency  u>  =  ub/Uco 
=  wing  first  torsional  mode  natural  frequency 
=  matrix  with  structural  frequency  ratios  squared 
along  main  diagonal 
(e.g.  (uji/uja)2,  (wm/wq)2) 

=  rrr11  modal  structural  coordinate  and  vector 
of  modal  structural  coordinates 


Subscripts 

0, 1  =  zeroth  and  first  harmonic,  respectively 

/  =  flutter  onset  (i.e.  neutral  stability)  condition 
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Chapter  8 

3D  Wing  With  Store:  F-16 
Wing/Store 

u  uiiimai  y 

•  Significant  progress  has  been  made  in  extending  out  reduced  order  method  (ROM)  and 

harmonic  balance  (HB)  method  for  computing  transonic  flutter  and  limit  cycle  oscilla- 
tions(LCO).  These  novel  solution  methods  provide  state  of  the  art  accuracy  using  either 
an  Euler  or  Navier- Stokes  computational  fluid  dynamics  (CFD)  fhodel  while  dramati¬ 
cally  decreasing  the  required  computational  time  and  effort.  This  new  capacity  has  been 
used  to  predict  the  effects  of  aerodynamics  nonlinearities  on  transonic  flutter  and  LCO 
of  a  F-16  wing/store  system  using  the  inviscid  (Euler)  and  viscous  (RANS)  versions 
of  the  HB  method.  * 

•  With  reference  to  the  attached  results  in  the  following  figures: 

Figure  8.1  shows  the  geometry  of  the  wing  and  a  portion  of  the  CFD  computational 
grid.  Mach  number  contours  for  the  steady  flow  about  the  wing  are  shown  in  Figure  8.2. 
The  reduced  velocity  at  which  flutter  occurs  and  the  eigenmodal  participation  factors 
at  the  flutter  condition  are  shown  in  Figure  8.3.  Finally,  in  Figure  8.4,  the  LCO 
respionse  for  the  wing  tip  acceleration  measured  in  g  is  shown  for  both  the  inviscid  and 
viscous  flow  models.  Note  that  although  the  reduced  velocity  at  which  flutter  begins  is 
similar  for  either  flow  model,  the  LCO  response  is  quite  different  depending  on  which 
flow  model  is  used,  inviscid  or  visous.  The  results  from  the  viscous  model  are  in  better 
agreement  with  the  LCO  response  observed  in  flight  tests. 

These  calculations  are  preliminary  quantitative  comparisons  between  the  theoretical 
results  and  flight  test  results  show  good  agreement  for  both  the  reduced  velocity  at  flutter 
and  the  flutter  frequency  as  well  as  the  general  level  of  LCO  response. 

•  More  detail  is  presented  in  the  following  section  (an  extended  abstract  for  AIAA/SDM 
2004).  For  additional  details  on  the  Navier-Stokes  (RANS)  calculations  and  results 
for  the  F-16  wing/store  LCO,  one  is  referred  to  the  paper  under  the  same  title  at  the 
upcoming  AIAA/SDM  2004  meeting  in  Palm  Springs. 

•  Presented  is  an  investigation  of  modeling  of  limit  cycle  oscillation  behavior  of  the  F-16 
fighter  configuration  using  a  frequency- domain  harmonic-balance  based  computational 
fluid  dynamic  approach.  Strong  nonlinear  limit  cycle  behavior  is  observed  in  the  case 


117 


of  a  viscous  flow  model  as  compared  to  an  inviscid  flow  model.  Details  of  the  compu¬ 
tational  model  and  plans  for  the  final  paper  are  presented. 


8.1  Introduction 

The  following  article  reports  the  latest  research  efforts  in  developing  and  demonstrating  the 
modeling  of  limit  cycle  oscillation  (LCO)  behavior  using  a  frequency  domain  harmonic  bal¬ 
ance  technique,  in  this  case  as  applied  to  a  F-16  fighter  aircraft  configuration. 

Denegri  [87]  has  presented  experimental  data  for  three  different  weapon  and  store  configu¬ 
rations  of  the  F-16  aircraft.  The  LCO  behavior  was  found  to  be  unique  for  each  of  these 
configurations. 

A  primary  purpose  of  this  investigation  is  to  study  the  possibility  of  the  using  the  harmonic 
balance  LCO  solution  method  to  model  some  of  the  various  LCO  trends  that  Denegri  [87] 
observed  for  the  limit  cycle  oscillations  of  the  F-16  wing. 

The  aerodynamics  of  the  F-16  wing  (but  not  the  stores)  are  modeled  by  either  an  Euler 
or  Navier-Stokes  based  CFD  code  and  the  solutions  are  obtained  using  the  novel  Harmonic 
Balance  Method  developed  at  Duke  University.  The  structural  model  contains  the  mass  and 
stiffness  characteristics  of  both  the  wing  and  the  stores. 


8.2  Theoretical  Details 

As  demonstrated  in  Thomas  et.  al.  [88,  89],  the  governing  aeroelastic  system  of  equations 
for  the  F-16  wing  can  be  written  in  the  frequency  domain  as: 


+  ifi)«  -  ce((,a)  =  o. 


(8.1) 


In  this  instance,  We  consider  a  linear  structural  model  expressed  in  terms  of  structural  nat¬ 
ural  modes  together  with  a  nonlinear  aerodynamic  model.  Structural  damping  is  neglected. 

kw  is  a  constant  dependent  on  the  wing  shape  and  overall  airplane  mass  given  by 

kw  _  7rAR(l+At)(l+At+At2) 

w  6m  ’  \  •  ) 

and  Cq  is  the  vector  of  generalized  aerodynamic  forces,  the  mfi1  element  of  which  is  given  by 

cQrn  =  — U  /  f  Pi  m  •  ndA  (8.3) 

Hoo  cr  J  JA 

where  the  integral  is  evaluated  over  the  surface  of  the  wing. 
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In  order  to  solve  the  aeroelastic  system  (Eq.  8.1),  one  must  be  able  to  compute  the  general¬ 
ized  aerodynamic  forces  C  q  for  finite  values  of  the  Structural  modal  coordinates  £.  This  is 
done  using  the  harmonic  balance  solver.  Details  of  the  harmonic  balance  procedure  can  be 
found  in  Hall  et.  al.  [90]  and  Thomas  et.  al.  [91,  88] 

8.3  New  Theoretical  Details 

M(-u;2I  +  fi2)£-Q.  (8.4) 

M(-o;2I  +  fl2)£-  QCMoo^u),^)  =  0.  (8.5) 

where 

<2i 

Q2 
Qm 


and 


or 

Qm==  'fidA 

(8.7) 

Qm  =  JJj>l  4>m  ■  ndA 

(8.8) 

Qm  =  Jj^  [pi  ($m  ‘  d  Ao)  +  P°  ($m  '  d  Ai)] 

(8.9) 

8.4 

LCO  Solution  Procedure 

As  far  as  the  LCO  solution  procedure,  we  begin  by  defining  the  residual  of  the  aeroelastic 
system  of  equations  (8.5)  as 

R(L)  =  +  n2)4/ft  -  2/6  =  0 

(8.10) 

where 
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Figure  8.1:  F-16  Inviscid  Flow  Case  Computational  Grid. 


L  =  \ 


Moo 

u 

Re(&)/6 


Re(W/6 


(8.11) 


8.5  Computational  Details 

For  the  results  presented  in  this  study,  the  CFD  method  is  a  variant  of  standard  Lax-Wendroff 
scheme  [92,  93,  94]  in  conjunction  with  the  one-equation  turbulence  model  of  Spalart  and 
Allmaras  [95]. 

Shown  in  Fig.  8.1  is  the  inviscid  computational  grid  near  the  F-16  wing  surface  and  wing 
symmetry  plane.  The  mesh  uses  65  (mesh  i  coordinate)  computational  nodes  about  the  wing 
in  the  stream-wise  direction,  33  (mesh  j  coordinate)  nodes  normal  to  the  wing,  13  (mesh  k 
coordinate)  nodes  along  the  semi-span,  and  12  more  nodes  wrapping  around  the  wing  tip. 
The  outer  boundary  of  the  grid  extends  five  wing  semi-spans  from  the  mid-chord  of  the  wing 
at  the  symmetry  plane.  A  similar  grid  is  used  for  the  viscous  flow  model  with  the  difference 
being  that  the  viscous  mesh  has  mesh  points  highly  clustered  near  the  wing  boundary. 
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Figure  8.2:  F-16  Inviscid  Flow  Case  Computed  Mach  Contours,  =  0.9,  a0  =  1.5  deg.. 


Figure  8.2  shows  computed  inviscid  Mach  number  contours  for  a  freestream  Mach  number 
of  Mqo  =  0.9  and  a  steady  angle-of-attack  of  o0  —  1*5  degrees. 

-X. 

In  the  following,  we  are  using  the  anti-symmetric  structural  mode  shapes.  For  the  aerody¬ 
namics,  we  are  modeling  only  one-half  of  the  wing.  %  As  such,  we  are  using  symmetric  unsteady 
aerodynamic  boundary  conditions.  So  as  an  initial  approximation,  we  are  assuming  that  the 
unsteady  aerodynamic  effects  from  one  side  of  the  overall  wing,  are  small  on  the  other  side 
of  the  wing.  For  the  final  paper,  our  plans  are  to  model  the  entire  wing  in  order  to  properly 
assess  the  magnitude  of  this  effect. 


8.5.1  F-16  Wing  Structural  Model 

Figure  8.3  shows  the  computational  nodes  of  the  structural  model  used  in  the  following 
analysis. 

The  structural  mode  shapes  correspond  to  the  “Typical  LCO”  configuration  of  the  F-16 
wing.  The  “Typical  LCO”  modes  shapes  are  for  the  case  of  the  F-16  with  the  following 
wing  stores:  Station  1:  Empty,  Station  2:  AIM-9L,  Station  3:  Air-Surface  Missile,  Station 
4:  Wing  Tank  (Empty).  See  Denegri  [87]  for  further  details. 

As  mentioned  previously,  only  the  anti-symmetric  mode  shapes  have  been  used  in  the 
HB/LCO  method.  For  this  initial  analysis,  the  following  results  are  based  on  using  the 
first  two  anti-symmetric  structural  modes.  • 
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Longitudinal  Coordinate,  x  (inches) 

Figure  8.3:  F-16  Structural  Model  Computational  Nodes 

8.5.2  Limit  Cycle  Oscillations 

Figure  8.4  shows  the  computed  viscous  flow  model  LCO  wing  tip  leading  edge  acceleration 
response  measured  in  g’s  for  three  difference  flight  altitudes.  The  current  model  does  not 
predict  any  significant  stable  or  unstable  LCO  behavior.  The  nearly  vertical  line  LCO  re¬ 
sponse  indicates  that  as  soon  as  flutter  onset  Mach  number  is  reached,  any  further  increase  in 
reduced  velocity  will  result  immediately  in  large  amplitude  unsteady  motions.  The  authors 
believe  there  are  three  possible  reasons  for  the  lack  of  strong  LCO  behavior  being  predicted 
by  the  HB/LCO  model. 

First,  there  could  be  issues  related  sufficient  mesh  resolution.  The  authors  wish  to  stress 
that  the  results  shown  are  preliminary.  Mesh  refinement  studies  are  currently  underway. 
Previous  work  by  the  authors  has  shown  that  sufficient  mesh  resolution  is  an  important 
aspect  of  accurate  LCO  prediction.  Second,  wing  external  stores  are  not  being  modeled  with 
the  current  CFD  method.  Other  researchers  have  suggested  that  unsteady  nonlinear  aero¬ 
dynamic  effects  brought  about  by  external  store  may  be  significant.  Third,  the  structural 
model  used  is  the  present  analysis  is  linear.  For  the  actual  F-16  configuration,  nonlinear 
structural  effects  could  play  a  significant  role  in  overall  LCO  behavior.  All  these  possibilities 
for  the  lack  of  strong  LCO  response  are  currently  being  investigated. 


8.6  Conclusions 

The  computation  of  limit  cycle  oscillation  behavior  of  a  F-16  fighter  aircraft  wing  configu¬ 
ration  is  demonstrated.  The  limit  cycles  are  computed  using  a  frequency  domain  harmonic 
balance  technique.  Calculations  are  ongoing.  Continuing  refinements  of  the  computational 
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Mach  Number, 

Figure  8.4:  Computed  F-16  LCO  Trends  for  Three  Difference  Altitudes  Based  on  Viscous 
Flow  Model.  a0  =  1.541  (deg) 


model  will  be  reported,  in  addition  to  further  details  of  the  overall  computational  method¬ 
ology.  ' 


Nomenclature 

Ar  =  aspect  ratio  =  wing  span  squared/wing  area 

b,c  =  wing  root  semi-chord  and  chord,  respectively 

C  q  =  vector  of  non-dimensional  generalized  forces 

I  =  identity  matrix 

fri  =  mass  of  wing 

Mqo  =  free- stream  Mach  number 

M  =  number  of  structural  modes 

At  =  generalized  mass  matrix 

pi  =  first  harmonic  unsteady  pressure 

Q  =  vector  of  generalized  aerodynamic  forces 

<7oo  =  free-stream  dynamic  pressure 

Re oo  =  free-stream  Reynolds  number 

Poo  =  free-stream  density 

Ua o  =  free-stream  velocity 

v  =  volume  of  a  truncated  cone  having  stream-wise 
root  chord  as  lower  base  diameter,  stream-wise 
tip  chord  as  upper  base  diameter,  and  wing 
half  span  as  height 

V  =  reduced  velocity  V  =  Uoo/ y/jhoab 

oo  =  wing  steady  (or  mean)  flow  angle-of-attack 
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At  =  wing  taper  ratio  \t  =  Cf/Cj. 
p  =  mass  ratio,  p  =  m/pooV 
ipm  —  rrfi1  structural  mode  shape 
a J,u>  =  frequency  and  reduced  frequency  u>  —  ujb/U^ 
Q  =  matrix  with  structural  natural  frequencies 
along  main  diagonal  (e.g.  o>i,  ui2,  ...,  u>m) 

£m,  £  =  modal  structural  coordinate  and  vector 
of  modal  structural  coordinates 

Subscripts 
r  =  root 

0, 1  =  zeroth  and  first  harmonic,  respectively 
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Chapter  9 

Commercialization  Strategy:  Phase  II 


With  a  strong  technical  background,  ZONA  also  has  extensive  experience  in  commercial¬ 
ization  of  its  product.  In  fact,  ZONA  has  been  serving  the  aerospace  community  through 
consulting/contractual  work  and  software  licensing  support  since  1986. 

ZONA’s  unsteady/steady  aerodynamic  software  in  the  aerospace  market  most  notably  in¬ 
cludes  ZONA51  (currently  the  Aero  Option  II  in  MSC/NASTRAN  -  an  industry  standard) 
and  the  ZAERO  software  system  (covering  the  entire  Mach  range).  ZONA  codes  have,  thus 
far,  accumulated  over  120  users  worldwide. 

Under  an  AF/STTR  Phase  II  contract,  ZONA  and  MSC  Software  are  reaching  a  busi¬ 
ness  agreement  to  commercialize  and  jointly  market  ASTROS*  (the  seamlessly  integrated 
ZAERO  module  into  ASTROS).  ASTROS  is  a  popular  Automated  STRuctural  Optimization 
System  software  among  aerospace  industry  and  universities. 


9.1  Next  Generation  Aeroelastic  Software 

After  some  20  years  of  continuing  R&D  effort  in  unsteady  transonic  aerodynamics,  high- 
level  (3-D  Euler  and  Navier-Stokes)  unsteady  CFD  methods  remain  inadequate  to  handle 
the  fixed-wing  transonic  aeroelastic  problems  encountered  in  industry.  ZONA’s  recent  mar¬ 
keting  survey  reveals  the  possible  causes  of  such  an  inadequacy.  These  include  that:  (1) 
the  current  CFD  methods  is  computationally  very  expensive  for  flutter  calculations,  (2)  the 
procedures  are  not  easily  user-adaptive,  (3)  Structure  engineers  prefer  to  work  with  solutions 
in  frequency-domain  not  time-domain.  Industry  standard  methods  such  as  the  classical  AIC 
method  is  lacking  in  the  transonic  regime  where  a  reliable  high-level  CFD  method  such  as 
CFL3D  is  badly  needed.  The  proposed  frequency-domain  POD/ROM  method  originated 
by  Hall/Dowell  has  a  great  potential  to  be  such  transonic  AlC-like  method  which  has  the 
essence  of  efficiency,  accuracy,  modularity  and  could  assess  the  key  physics  for  nonlinear 
flutter/LCO.  Its  future  applicability  is  far  reaching  which  extends  to  areas  in  aeroelasticity, 
aeroacoustics,  turbomachinery,  ASE,  MDO,  etc.  On  the  other  hand,  the  ZAERO  aeroelastic 
system  exclusively  adopts  the  AIC  methodology  as  basis  for  the  unified  unsteady  aerody¬ 
namics.  Although  efficient,  ZAERO  is  incapable  to  provide  nonlinear,  viscous  unsteady 
aerodynamics  with  large  amplitude  oscillations.  Once  developed,  the  transonic  frequency- 
domain  POD/ROM/HB  (Harmonic  Balance)  methodology  along  with  its  Navier-Stokes  op¬ 
tion,  would  be  an  integral  part  of  ZONA’s  ZAERO  aeroelastic  system.  ZONA  envisions  that 
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the  production  code  of  the  frequency-domain  POD/ROM/HB  method  should  become  the 
future  industry  standard  for  routine  transonic  flutter/LCO  analysis. 


9.2  Commercialization  of  the  Software  with  ZAERO 

As  the  frequency-domain  POD/ROM/HB  code  is  successfully  developed,  ZONA  intends  to 
package  the  POD/ROM/HB/Navier-Stokes  methodology  (with  either  software  option)  into 
a  well-defined  commercial  software  product.  To  facilitate  its  sales,  ZONA  plans  to  couple  the 
code  with  the  unified  aerodynamic  module  within  ZAERO  and  jointly  market  for  both  codes. 

ZONA  has  currently  teamed  up  with  a  major  FEM  software  house  (MSC  Software)  for 
ZONA51  (in  MSC/NASTRAN)  and  ZAERO  licensing.  ZONA  currently  supports  over 
120  users  worldwide  for  the  MSC/NASTRAN  Aero  Option  II.  Released  in  April  1999, 
ZAERO  has  already  had  a  dozen  users  including  Lockheed  Martin/Vought  Systems,  Boe¬ 
ing/Commercial,  DLR/Gttingen,  SDRC,  NASA/MSFC,  NASA/Dryden,  Scaled  Composites 
Coleman  Aerospace,  etc.  With  ZONA’s  extensive  aerospace  contractual/software  licensing 
experience,  this  software  product  is  expected  to  break  into  the  aerospace  market  quite  read¬ 
ily.  Prospective  customers  include  ZONA’s  current  clientele,  other  aerospace  companies  and 
DoD  organization. 

As  a  first  license  precurser  /3-sites  will  be  set  up  with  Boeing/St.  Louis  and  Lockheed 
Martin/Forth  Worth  to  upgrade  the  code  as  a  viable  industry  production  software. 

Other  commercialization  steps  include  advertising  through  magazines,  the  internet,  and  at¬ 
tending  AIAA  conventions.  All  these  are  currently  being  pursued  by  ZONA  for  the  licensing 
sales  of  ZONA  codes. 
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Chapter  10 
Conclusion 


ZONA  Technology,  Inc.  and  Duke  University,  the  ZONA/Duke  team,  has  completed  a 
successful  Phase  II  program  entitled  ’’Nonlinear  Reduced-Order  Modeling  of  Limit  Cycle 
Oscillations  of  Aircraft  Wings  and  Wing/Stores” .  The  overarching  goal  for  the  Phase  II 
work  is  to  build  on  the  very  substantial  progress  made  in  Phase  I  and  to  develop  a  capa¬ 
bility  for  a  highly  computationally  efficient  and  physically  accurate  mathematical  modeling 
of  limit  cycle  oscillations  (LCO)  and  other  nonlinear  aeroelastic  phenomena.  The  approach 
uses  the  concepts  of  aerodynamic  modes  as  well  as  structural  modes.  Such  models  provide 
substantially  improved  physical  understanding  and  also  more  accurate  prediction  of  LCO  in 
particular  and  nonlinear  aeroelastic  responses  in  general.  This  capability  may  also  lead  to  a 
rational  prediction  of  buffet  onset  due  to  aerodynamically  nonlinear  oscillations. 

The  specific  objective  of  the  proposed  program  is  to  develop  the  frequency-domain  POD/ROM/ 
HB  EigenMode  method  (HB  method)  and  to  generalize  this  method  to  three-dimensional 
viscous  transonic  flow,  in  order  to  further  the  understanding,  provide  accurate  and  efficient 
prediction,  and  exercise  control  of  LCO/flutter  for  aircraft  wings  and  wing/stores.  In  so 
doing,  the  high-level  Reynolds-average  Navier-Stokes  (RANS)  Equation  is  solved  by  the  de¬ 
veloped  HB  method  framework  for  computational  expediency,  while  maintaining  solution 
accuracy  for  LCO/flutter  prediction. 

The  HB  method  developed  is  aimed  at  the  capability  of  handling  the  following  physical 
effects: 

•  Three  dimensional  (3D)  as  well  as  2D  flows. 

•  Nonlinear  as  well  as  time  linearized  aerodynamic  pressures  and  forces. 

•  Multi-body,  e.g.,  wings  plus  stores,  as  well  as  single  body  (e.g.,  wing  only)  configura¬ 
tions. 

•  Many  as  well  as  few  structural  modes. 

•  LCO  (nonlinear)  as  well  as  flutter  (linear)  aeroelastic  prediction. 

Accordingly,  the  ZONA/Duke  team  has  successfully  accomplished  the  following  technical 
objectives: 

•  Develop  the  2-D  N-S  version  of  the  above  proposed  approach. 
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•  Investigate  further  the  intial  condition  impart  on  LCO  via  a  2-D  time-domain  method 
(CFL3D/N-S  code). 

•  Extend  the  geometry  capability  to  3-D  including  wing/store  configurations. 

•  Generalize  the  2-D/3-D  methods  to  include  Harmonic  Balance  technique  for  nonlinear 
LCO  analysis. 

•  Include  multiple  structural  mode  analysis  capability. 

•  Perform  LCO  case  studies  including  wing-alone  and  F-16  wing/store  systems. 

Throughout  the  course  of  the  HB  method  development,  the  following  conclusions  from  the 
tasks  (presented  in  Chapter  2  through  Chapter  8)  are  worth  noting: 

1.  Transonic  LCO/flutter  solutions  using  the  2D  HB  method  for  inviscid  and  viscous 
flows:  Frequency-Domain. 

•  An  extensive  parameter  analyses  of  an  airfoil  aeroelastic  configuration  (NACA64A010A 
airfoil)  has  been  achieved  using  a  highly  efficient  time  linearized  CFD  computa¬ 
tional  technique.  Correlation  with  experiment  remains  an  open  challenge,  how¬ 
ever  a  comparison  of  results  from  our  CFD  model  with  results  from  an  unsteady 
aerodynamic  experiment  by  Davis  and  Malcom[37]  for  lift  and  moment  shows 
encouraging  agreement. 

•  An  efficient  procedure  for  computing  the  LCO  behavior  of  aeroelastic  airfoil  con¬ 
figurations  (NLR7301  airfoil)  in  viscous  transonic  flows  is  presented.  Viscous 
effects  are  demonstrated  to  be  an  important  factor  in  determining  the  behavior  of 
LCO  response  with  respect  to  reduced  velocity.  Viscous  effects  for  the  NLR  7301 
aeroelastic  configuration  under  investigation  lead  to  a  much  more  gradual  rate 
of  increase  in  LCO  amplitude  with  respect  to  reduced  velocity  when  compared 
to  an  inviscid  model.  The  rapid  increase  in  the  computed  LCO  response  beyond 
the  flutter  point  emphasizes  the  importance  of  making  careful  and  comprehensive 
experimental  measurements  of  LCO  response  over  a  range  of  reduced  velocities 
or  other  parameters  as  the  flutter  boundary  is  exceeded. 

2.  Transonic  LCO/flutter  solutions  using  CFL3D.AE  code:  Time-Domain 

•  A  CFD  time-marching  method,  CFL3D  v6,  is  employed  to  investigate  the  in¬ 
fluences  of  viscosity,  initial  airfoil  position,  and  initial  perturbation  on  the  am¬ 
plitudes  of  the  transonic  limit  cycle  oscillations  (LCO)  of  a  supercritical  airfoil. 

As  expected,  stronger  aerodynamic  nonlinearity  leads  to  smaller  LCO  amplitudes, 
even  a  damped  solution  while  weaker  aerodynamic  nonlinearity  incurs  larger  LCO 
amplitudes,  even  a  divergent  solution. 

•  The  computed  results  due  to  different  methods  on  the  transonic  LCO  solution  for 
NLR7301  supercritical  airfoil  at  M  =  0.75  can  be  summarized  by  the  following 
figure  prepared  by  Duke  University.  It  is  seen  that  the  DLR  measured  LCO 
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amplitude  (Schewe)  is  only  a  fraction  of  the  computed  results  including  that  of 
ZONA  (Tang  et  al),  Naval  Postgraduate  School  (Weber  et  al)  and  Duke  University 
(Hall  and  Dowell,  Hollow  Symbols). 

•  It  is  clear  that  the  HB  method  is  superior  to  the  time-domain  methods  such  as 
CFL3D  in  terms  of  computational  efficiency  and  the  ease  of  LCO/flutter  solution 
identification. 

3.  Transonic  LCO/flutter  solutions  for  3D  wings  and  wing/store  systems  using  the  HB 
method:  Frequency-Domain 

•  A  harmonic  balance  method  for  modeling  nonlinear  periodic  unsteady  three- 

dimensional  inviscid  transonic  flows  about  wing  configurations  (AGARD  445.6 
wing)  is  presented.  Demonstrated  is  the  ability  of  the  method  to  model  nonlinear 
aerodynamic  effects  due  to  large  amplitude  motions  for  a  well  known  aeroelastic 
configuration.  The  use  of  only  two  harmonics  in  the  procedure  is  shown  to  be 
quite  sufficient  for  producing  harmonic  convergence.  A  limit  cycle  oscillation  so¬ 
lution  methodology  is  also  presented  and  demonstrated  for  the  same  benchmark 
transonic  configuration.  A  weak  LCO  is  observed,  which  is  consistent  with  results 
from  time-domain  calculations  by  Gordnier  and  Melville  [86].  Nevertheless,  it  is 
suggested  that  nonlinear  aerodynamic  effects  may  at  least  partially  explain  some 
of  the  previously  observed  differences  between  experiment  and  linear  aeroelastic 
theory.  •••* 

•  Significant  progress  has  been  made  in  extending  out  reduced  order  method  (ROM) 
and  harmonic  balance  (HB)  method  for  computing  transonic  flutter  and  limit 
cycle  oscillations  (LCO).  These  novel  solution  methods  provide  state  of  the  art 
accuracy  using  either  an  Euler  or  Navier-Stokes  computational  fluid  dynamics 
(CFD)  model  while  dramatically  decreasing  the  required  computational  time  and 
effort.  This  new  capacity  has  been  used  to  predict  the  effects  of  aerodynamics 
nonlinearities  on  transonic  flutter  and  LCO  of  a  F-16  wing/store  system  using  the 
inviscid  (Euler)  and  viscous  (RANS)  versions  of  the  HB  method. 

Finally,  this  Phase  II  final  report  intends  to  provide  a  comprehensive  description  of  our 
approach  and  results  obtained  in  achieving  the  Phase  II  technical  objectives.  These  achieve¬ 
ments  provide  a  sound  foundation  upon  which  the  Phase  III  effort  for  technology  transfer 
for  a  production-ready  LCO  prediction  software  can  be  carried  out  thus  rendering  viable 
commercialization. 
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